AD-A064  676  AEROSPACE  CORP  EL  SEGUNOO  CA  ENGINEERING  GROUP 

THE  ECLECTIC  SIMULATOR  PROGRAM  (ESP)  USAGE  GUIOE.(U) 

MAT  80  E  R  COFFEY »  H  J  WERTZ  F04701-79-C- 

UNCLASSIFIED  TR-0080(9320)-l  SD  -TR-80-21 


1  c*3 


REPORl  SD-TR-80-21 


CO 

<>• 

CO 

oc 

o 

c 


Usage  Guide 


EMMAGENE  R.  COFFEY 
in  Association  with 
HARVEY  J.  WERTZ 
Mission  Information  Systems  Division 
Engineering  Group*  f 
The  Aerospace  Corporation 
El  Segundo,  Calif.  90245 


1  May  1980 


Interim  Report 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED 


CJ> 


UJ 


THE  AEROSPACE  CORPORATION 

Prepared  for 

SPACE  DIVISION 
AIR  FORCE  SYSTEMS  COMMAND 
LOS  ANGELES  AIR  FORCE  STATION 
P.O.  Box  92960,  Worldway  Postal  Center 
Los  Angeles,  Calif.  90009 


MO  5  27  078 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Date  Entered) 


(J  X?  REPORT  DOCUMENTATION  PAGE 


■  REPOST  NUMBER 

s  dVt  R  -  8/b  -  2  if 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


A. 


[2.  GOVT  ACCESSION  NO. 


J-  RECIPIENT'S  CATALOG  NUMBER 


tfffflj'"'  * 


[  (*)  JTHE  ECLECTIC  SIMULATOR  PROGRAM  (ESP)/ 
Lfl  USAGE  GUIDE <  *  . .  .  -  . .  ‘ —  * 

- . 


S.  TYPE  OF  REPORT  A  PERIOD  COVERED 

Interim  Report 
7-1-75  to  12-1-79 


4 


T  R-00  80(9  320 )-!/ 


NUMBER 


7.  AUTHORf.; 


’  Eramagene  R. ^Coffey, ^Harvey  J.^Wertz 

10 


S.  CONTRACT  OR  GRANT  NUMBE 


1 -79- C- 


BEIVJJ 

f 


RFORMING  ORGANIZATION  NAME  AND  ADDRESS 

The  Aerospace  Corporation 
El  Segundo,  California  90245 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  S  WORK  UNIT  NUMBERS 


({&  W 


11.  CONTROLLING  OFFICE  NAME  AND  ADORESS 


12-  REPORT  DATE - 

< i iy.  i Myf 1 i°  / 

Vv__^TTS.  Number  OF  PAGES 
186 


U  MONITORING  AGENCY  NAME  A  ADDRESS^!  dllltrtnl  from  Controlling  Olllcm) 

Space  Division 

Air  Force  Systems  Command 
Los  Angeles,  Calif.  90009 


IS.  SECURITY  CLASS,  (of  thlo  report; 

Unclassified 


i"s».  DECLASSIFl cation/ downgrading 
schedule 


16.  DISTRIBUTION  STATEMENT  fof  ffiio  Roporf; 

Approved  for  public  release;  distribution  unlimited. 


17,  DISTRIBUTION  STATEMENT  (of  lire  mbttrmct  entered  In  Block  30.  II  different  Horn  Report; 

7$  -X: '  J  • 


A  ' 


/  — ..  '  r  -  J 


y 
■  J 


8.  SUPPLEMENTARY  notes 


19.  KEY  WOROS  (Continue  on  re  veree  at  do  1/  neceeeery  end  Identity  by  block  ntmtbmr) 


Simulation 
Integration 

Numerical  Integration 
Ordinary  Differential  Equations 
D£scontinuous  Driving  Functions 


Nu 


Hysteresis  Nonlinearities 
Computer  Simulation 
Predictor-Corrector  Integration 
Runge-Kutta  Integration 
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The  Eclectic  Simulator  Program  (ESP)  is  a  system  (a  precompiler  plus  a 
collection  of  subroutines)  that  permits  the  fast,  easy  solution  of  ordinary 
differential  equations.  Any  user  with  a  general  knowledge  of  F0RTRAN  can 
utilize  ESfs  many  labor-saving  devices  to  code  a  problem  with  minimal 
effort.  Special  ESP  features  permit  translation  of  engineering  blocks,  dis¬ 
continuities,  and  hysteresis  patterns  directly  into  computer  code,  and  the  use 
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of  WHELP  in  conjunction  with  ESP  facilitates  efficient  coding  of  matrix 
algebra  equations.  Simple  input  cards  enable  the  user  to  directly  control 
solution  and  timing  accuracy  and  to  specify  or  change  run  times,  initial 
conditions,  and  various  other  parameters  easily  when  making  multiple  or 
stacked  runs.  Finally,  ESP  allows  the  user  to  select  from  a  wide  variety 
of  output  options.  This  manual  is  intended  to  be  both  a  learning  tool  for 
the  novice  and  a  detailed  reference  for  the  experienced  user.  i 
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PREFACE 


This  revision  of  The  Eclectic  Simulator  Program  Usage  Guide  was 
prompted  by  a  number  of  modifications  and  improvements  to  the  ESP  package 
which  have  been  made  since  1975.  The  primary  reason  for  the  changes  was 
the  desire  to  have  common  versions  of  ESP  and  WHELP  available  on  both  the 
CDC  7600  and  the  IBM  370.  From  the  user's  point  of  view  this  has  been 
achieved  even  though  both  WHELP  and  PREC$MP  are  written  in  PLI  for  the 
IBM  370  and  in  F0RTRAN  for  the  CDC  7600,  There  were,  however,  a 
number  of  changes  required  which  will  affect  existing  programs.  The  first  of 
these  is  a  revision  of  variable  and  common  block  names.  The  second  is  a 
revision  and  restructuring  of  the  major  subroutines  to  simplify  program  flow 
and  logic  and  to  improve  error  control  and  handling  of  potential  job  abort 
situations.  The  third  is  replacement  of  the  former  default  integration  algo¬ 
rithm  with  the  Shampine  variable -order /variable  step  routine,  which  appears 
to  outperform  the  Hamming  Predictor  Corrector  in  both  speed  and  accuracy 
on  a  variety  of  applications  tested  thus  far. 

In  effecting  these  changes  to  the  ESP  package,  considerable  effort  has 
been  made  to  minimize  changes  required  of  the  user.  The  only  pervasive 
changes  affecting  the  user  are  the  shortening  of  variable,  subroutine,  and 
common  block  names  to  six  characters  or  less  and  some  rearranging  of 
common  blocks,  both  necessitated  by  IBM  FORTRAN  constraints.  In  general, 
the  contents  of  this  guide  reflect  a  CDC  orientation  with  IBM  counterparts 
noted  where  possible.  However,  an  ESP  program  written  according  to  this 
guide  may  be  run  interchangeably  on  either  IBM  or  CDC  by  simply  changing 
the  appropriate  control  cards. 

In  addition  to  revisions  of  all  portions  of  this  guide  reflecting  the  above 
changes,  the  appendix  relating  to  WHELP  (Appendix  H)  has  been  greatly 
expanded  to  include  all  current  capabilities  of  WHELP.  Notice  that  WHELP 
may  be  used  with  or  without  ESP,  and  thus  Appendix  H  may  be  used  without 
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reference  to  the  remainder  of  the  guide  and  is  in  fact  the  most  complete 
current  documentation  of  WHELP. 

Hopefully,  this  will  be  the  final  major  rewrite  of  both  the  software 
ana  the  documentation.  In  particular,  this  edition  of  the  manual  has  been 
designed  for  easy  updating  when  the  inevitable  (wishfully  small)  changes  are 
made  to  the  software. 
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SECTION  I 


INTRODUCTION 


The  Eclectic  Simulator  Program  (ESP)  enables  the  user  to  solve 
ordinary  differential  equations  with  speed,  accuracy,  versatility,  and  mini¬ 
mal  effort.  The  user  codes  only  that  information  unique  to  his  particular 
problem:  the  differential  equations,  initial  conditions,  and  desired  output. 
This  information  must  be  coded  in  ESP  language,  which  is  a  special  purpose 
programming  shorthand  developed  just  for  this  program,  and  documented  in 
this  manual.  ESP  then  does  literally  all  the  remaining  work. 

To  accomplish  its  purpose,  the  ESP  system  is  composed  of  two  parts-- 
a  precompiler  and  a  set  of  F0RTRAN  subroutines.  The  precompiler  reads 
the  ESP  shorthand  code  written  by  the  user  (This  code  will  not  look  like  a 
F0RTRAN  program.  )  and  translates  it  into  F0RTRAN,  while  adding  the 
necessary  cards  (such  as  C0MM0N  block,  DIMENSION  statements,  and 
RETURNS)  to  produce  complete  and  executable  F0RTRAN  subroutines.  The 
output  of  this  precompiler  is  then  joined  with  the  second  part  of  the  ESP 
package,  the  subroutines  which  do  the  integration  and  other  chores,  to  forma 
complete  F0RTRAN  program  which  is  then  executed  by  the  computer. 

ESP  is  not  only  highly  efficient  in  terms  of  both  user's  effort  and  com¬ 
puter  running  time,  but  it  is  also  highly  flexible.  A  number  of  special  capa¬ 
bilities,  labor  saving  devices,  and  alternate  means  to  the  same  ends  are  part 
of  the  package;  however,  the  user  has  considerable  latitude  in  deciding  which 
features  he  will  use  and  how  large  or  how  small  a  problem  he  wishes  to  solve. 
Briefly  listed  below  are  some  of  the  distinctive  features  and  capabilities  of 
ESP: 

•  The  derivatives  are  normally  defined  as  first-order  differential 

equations,  but  may  instead  be  translated  directly  from  engineering 
block  diagrams  to  *BL0CK  cards  without  any  intervening  algebra. 
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•  The  basic  integration  algorithm  is  the  highly  efficient  SHAMPINE' 
method,  which  combines  variable  stepsize  with  variable  order 
integration  in  response  to  continuous  error  checks,  but  the  user 
may  opt  to  run  ESP  using  any  of  several  other  integration  algo¬ 
rithms,  namely,  second-  or  fourth-order  Runge-Kutta,  or 
Predictor-Corrector  (which  uses  Runge-Kutta  as  a  starter),  any 
of  which  may  be  run  with  either  a  fixed  or  variable  stepsize. 

•  Significant  sign  changes,  discontinuous  driving  functions,  hyster¬ 
esis  nonlinearities  and  the  like  may  all  be  accurately  and  easily 
coded  into  the  system  by  means  of  special  ESP  language  command 
cards. 

•  Output  from  an  ESP  program  may  take  many  forms,  such  as  auto¬ 
matically  formatted  print,  user-formatted  print,  calcomp  pen 
plots,  printer  plots,  microfilm  plots,  or  magnetic  tape  files. 

•  Since  inputs  such  as  initial  conditions,  run  times,  and  parameters 
can  be  easily  changed,  a  series  of  runs  or  a  set  of  "stacked"  runs 
can  be  made  with  a  minimum  of  effort. 

•  The  user  can  directly  control  the  degree  of  accuracy  required  for 
the  problem  solution  and  also  for  the  timing  of  discontinuities 
with  simple  input  cards. 

•  Vector-Matrix  expressions  can  be  used  in  their  natural  form  to 
compute  derivatives,  by  using  WHELP  along  with  ESP. 

This  manual  attempts  to  meet  the  needs  of  both  the  novice  and  the  ex¬ 
pert  user  of  ESP.  It  is  hoped  that  sufficient  explanation  and  examples  have 
been  given  in  Sections  II  through  VII  to  enable  the  uninitiated  to  write  a  suc¬ 
cessful  program.  On  the  other  hand,  considerable  detail  has  been  included 
throughout  to  aid  all  users  in  answering  their  own  questions  and  debugging 
their  own  programs. 

Section  II  includes  a  straightforward  example  of  ESP  usage,  from 
problem  definition  through  printed  and  plotted  output.  Careful  study  of  this 
example  and  its  annotation  should  give  the  user  a  helpful  overview  of  how 
ESP  works  and  how  the  various  parts  of  user-coding  relate  to  each  other. 

Referred  to  in  this  manual  as  ADAMS  because  it  is  an  Adams -like  method 
and  SHAMPINE  is  too  long  for  a  FORTRAN  name 


Following  this  example,  material  is  arranged  topically  by  sections,  one 
section  for  each  major  aspect  of  setting  up  an  ESP  program:  defining  the 
derivatives,  selecting  the  integration  method,  modeling  discontinuities, 
specifying  output,  and  defining  input.  Each  section  begins  with  an  overview  of 
the  capabilities  relating  to  the  section  topic,  and  then  discusses  each  in  detail, 
starting  with  the  simplest  and  most  basic  usage  and  progressing  to  more 
complicated  options  and  considerations  near  the  end. 

It  is  strongly  recommended  that  the  user's  first  attempt  at  coding  ESP 
involve  a  fairly  simple  problem  or  a  simplified  version  of  a  larger  problem, 
and  that  more  complex  aspects  of  ESP  be  added  only  after  the  basic  ones  are 
well  understood  and  seem  to  be  working  properly.  In  keeping  with  this 
approach,  it  is  suggested  that  the  user  study  the  example  problems  (Section  II) 
and  read  the  first  few  pages  of  each  of  Sections  III  through  VII  before  attempt¬ 
ing  to  code  his  first  problem.  Later  parts  of  Sections  III  through  VII  and  the 
Appendices  may  be  regarded  more  as  reference  material  and  used  only  as 
needed,  although  particular  attention  should  be  called  to  Appendix  D,  Program 
Control  and  Execution,  for  those  who  wish  to  understand  more  fully  how 
ESP  works.  The  Table  of  Contents  and  Index  should  facilitate  easy  location 
of  any  other  material  of  interest. 
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II.  SOME  SIMPLE 
EXAMPLE  PROBLEMS 


SECTION  II 


SOME  SIMPLE  EXAMPLE  PROBLEMS 


To  aid  the  user  in  obtaining  an  overview  of  how  ESP  works,  this  section 
consists  solely  of  two  example  problems.  The  first  is  a  very  simple  or 
minimum  problem  and  includes  the  statement  of  the  problem  and  the  coding 
required  to  solve  i'  using  ESP.  The  second  example  is  slightly  more  com¬ 
plex  and  is  fully  discussed  from  problem  definition  through  analysis,  coding, 
and  resulting  printout. 

A.  EXAMPLE  1 

Prob lem:  Integrate  the  following  differential  equations  from  t  =  0  to 
t  =  10.0  sec,  printing  t,  Y,  and  Y  every  0,  5  sec: 

Y  =  cos  G  Y  +  Bt 

Y  =  Y  +  sin  0  Y 


where 


all  initial  conditions  =  0. 

6  =  0.  4*t 
B  =  0. 142 

Coding: 


[Control  cards--see  APPENDIX  B] 

*  DERI  VS 

THETA  =  0.  4  *  T 

DY  (1)  =  COS  (THETA )  *  Y(l)  +  0.  142  *  T 
DY (2)  =  DY (1)  +  SIN(THETA)  *  Y(l) 
*ENDDERIVS 

*PRINT  TIME  =  T  $  Y  =  Y(l)  $  YDOT  =  Y(2)  $  $ 
*RUN  2  0.  0.5  10.0  $ 


\ 
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B.  EXAMPLE  2 

This  example  will  be  presented  in  detail  from  problem  definition  through 
to  the  output  of  a  completed  program  in  five  steps.  Step  1  is  the  statement  of  a 
problem  as  the  typical  user  might  define  it.  Step  2  shows  a  step-by-step 
analysis  of  this  problem  and  translation  of  its  characteristics  into  ESP  code, 
while  Step  3  illustrates  the  actual  arrangement  of  this  code.  Steps  4  and  5 
are  provided  by  ESP  and  show  the  F0RTRAN  output  of  the  precompiler  and  the 
actual  printed  and  plotted  output  requested  by  the  user. 


STEP  1:  STATEMENT  OF  THE  PROBLEM 
Integrate  the  above  system  from  T=0,  to  T=0.  5 
Inputs:  Gj  =  5.0 

g2  =  1.0 

a  =0.1 


b  =0.01 
£  =0.5 

U)  =  1.0E1 

All  initial  conditions  =  0. 
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Outputs:  Print  every  0.02  second  the  following  values  and  labels: 

Time  =  T 

Error  =  1.0  -  G^Y^ 

Output  = 

Plot:  Output  versus  Time 

Error  versus  Time 

STEP  2:  ANALYSIS  OF  THE  PROBLEM 

1.  The  number  of  integrations  required  (3),  the  run  interval  and  the  print 
interval  will  be  specified  on  the  *RUN  card 

*RUN  3  0.  0  0.02  0.  5  $  [See  Section  VII-A] 

2.  Input  constants  will  be  input  on  a  *PAR  card  and  equivalenced  to  their 
names  in  the  equations  so  that  later  they  may  be  easily  changed. 

PAR(l)  «  Gj  *  G1  =  5. 

PAR (2)  =  b  *  B  =  0.01 
PAR  (3 )  =  a  =  A  =  0.  1 
PAR<4)  *  (  =  ZETA  =0.5 

PAR (5)  a  W  =  0MEGA  =  1.  0E1 
PAR (6)  a  G2  a  G2  =  1.0 

by  using  *PAR  5.  0.01  0.1  0.5  1.E1  1.0  $  [See  Section  VII-C] 

represents  the  following  characteristics: 

input  =  1.  0  -  G^  *  Y^ 

If  input  ^  0,  output  =  0. 
output  =  If  input  >  0,  output  =  1.0 

To  detect  the  exact  point  at  which  the  value  of  input  changes  sign  and  to 
set  the  proper  output,  the  *SWTCH  feature  will  be  used: 

*SWTCH  1  1.0  $  0.  $  1. 0-PAR(6)*Y(3)  $ 

[See  Section  V-A] 


2-3 


G1 

s  +  b 

is  equivalent  to  =  Gj  *input-b*Y(l)  where 
input  =  SWCH1  which  is  the  value  resulting 
from  the  *SWTCH  statement. 


This  will  be  coded  as 

DY(1)  =  G1*SWCH1  -  B * Y ( 1 )  [See  Section  III-A] 


is  a  second-order  block  which  is  equivalent  to 
Y3  =  Y^s+aJ/fs2  +  2{(js  +CJ2)] 

This  could  be  solved  for  Y^  in  terms  of  its  auxiliary  function  Y^  and  coded  as 

DY (2)  =  A*Y(1)  -  0MEGA**2*Y(3) 

DY (3)  =  -2.  O*ZETA*0MEGA*Y(3)  +  Y(2)  +  Y(l) 

but  it  is  faster  and  easier  to  code  it  by  using  the  *BL0CK  input  feature 

*  BLOCK  2  1.0  A  2.  *ZETA*0MEGA  0MEGA**2  Y(3)  Y(2)  Y(l)  $ 

[See  Section  III-C] 

6.  Printing  and  storing  of  plot  data  will  be  done  by  using  the  *PRINT 
statement 

*PRINT  TIME=PL0T(1)  =  T  $  ERR0R=PL0T(2)  =  1.  0-PAR(6)*Y(3)  $ 

0UTPUT=PL0T(3)  =  Y(3)  $  $  [See  Section  VI-A- 1] 


s  +  a 


s2  +  2f,ws  +<*? 


7.  Printer  plots  with  all  default  features  will  be  generated  by  using 
*GR A PH  1  3 

0UTPUT  VERSUS  TIME  [See  Section  VI-B-2] 

*GRAPH  1  2 

ERR0R  VERSUS  TIME 
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8. 


A  title  will  be  assigned  to  ail  output  pages  by  using  the  *TITLE  card 


^TITLE  EXAMPLE  FOR  ESP  MANUAL 

[See  Section  VII-G-2] 


Beginning  ot  FORTRAN  subroutines  written  by  the 
precompiler  from  the  user's  coding 

This  mein  progrem  is  normelly  written  entirety  by  PRCC0MP. 
It  the  user  writes  his  own.  it  should  be  pieced  et 
the  beginning  of  the  deck. 
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FROGRArt  LENGTH  12B  10  have  been  coded  by  the  user  on  * SWTCH  cat 

SCH  LABELED  COMMON  LENGTH  541B  353  The  output  of  these  switches  is  computed  in 

47000B  SCM  USED  SUBROUTINE  DERIVS. 


•  dummy  or  do-nothing  routine,  written  entirety 
by  PREC0MP. 
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PROGRAM  LENGTH  16B  14 

SCM  LA3ELED  COMMON  LENGTH  541B  353 

47000B  SCM  USEO 


LOADER  VERSION  1.0  06/07/79  13.54.35.  PAGE 
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06/07/79  EXAMPLE  FOR  ESP  MANUAL  PAGE  NO 


OUTPUT  VERSUS  TIME 


1.200E-01  2 .400E-0I  3.600E-01  4.800E-01  6.000E-01 


TIME  =  .056  SEC.  time  required  for  plotting 
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SECTION  III 


DEFINING  DERIVATIVES 


The  most  important  segment  of  coding  which  the  user  must  provide  is 
that  which  defines  his  derivative  equations.  This  segment  will  be  translated 
by  the  precompiler  into  SUBROUTINE  DERIVS  which  is  then  called  as 
needed  by  the  ESP  integration  package  for  each  evaluation  of  the  derivatives. 
The  user,  therefore,  must  code  his  equations  in  a  manner  that  can  be  recog¬ 
nized  and  properly  interpreted  by  ESP.  There  are  several  ways  to  do  this, 
and  the  user  can  choose  the  one  or  more  ways  most  suitable  to  his  problem 
from  the  following  alternatives: 

•  The  derivatives  may  be  written  as  a  set  of  first-order  differential 
equations  in  terms  of  the  ESP  variables  Y(i),  DY(i),  and  T. 

•  The  derivatives  may  be  written  in  terms  of  the  user's  variables 
and  their  values  then  "moved"  into  the  ESP  array,  DY.  (This 
option  is  particularly  desirable- -frequently  necessary- -if 
WHELP  expressions  are  used  to  compute  the  derivatives.) 

•  The  derivatives  may  be  written  as  *BL0CK  statements,  which 
permit  direct  and  simple  translation  of  engineering  block  diagrams 
directly  into  code  that  ESP  can  interpret. 

In  any  case,  all  derivatives  must  be  defined  within  a  section  of  coding 
which  begins  with  the  card  *DERIVS  starting  in  column  1  and  is  terminated 
with  the  card  *END DERIVS,  also  starting  in  column  1.  This  section  of  coding 
may  be  placed  first,  after  any  user  supplied  subroutines,  or  it  may  be  placed 
after  the  *ICC0MP.  . .  *ENDIC  or  *(&UTPUT. . .  *ENDtf)UT  sections.  (See 
Appendix  A-3  User's  Deck  Structure.  )  In  addition  to  the  derivatives,  this 
section  will  contain  any  *SWTCH  or  *SWMEM  cards  used  (See  Section  IV 
Discontinuities.  ). 
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A.  DEFINING  THE  DERIVATIVES  AS  FIRST-ORDER 

DIFFERENTIAL  EQUATIONS 

SUBR0UTINE  DERIVS  is  written  by  ESP  from  the  segment  of  the  user's 

coding  beginning  with  the  :;:DERIVS  card  and  ending  with  the  *ENDDERIVS  card. 

It  will  receive  from  the  integration  package  the  value  of  the  independent 

th 

variable  T  and  the  vector  Y  containing  the  values  of  each  i  integral  and  must 
return  to  it  the  vector  DY  containing  the  derivatives  of  each  corresponding 
Y(i).  Therefore,  certain  points  must  be  kept  in  mind  in  coding  the  derivative 
equations : 

•  Each  equation  must  be  solved  for  some  Y  so  that  it  can  be 
written  in  the  form  DY(i)  =  some  expression,  one  DY(i)  per 
integrator  required. 

•  The  variable  T  is  always  the  independent  variable.  (Its  range  of 
values  is  specified  on  the  *RUN  card  explained  in  Section  VII- A.) 

•  Y(i)  should  always  be  assumed  to  represent  the  result  of  integra¬ 
tion  at  the  current  value  of  T. 

•  DY(i)  may  appear  on  the  right-hand  side  of  an  equation  if  it  has 
been  defined  above,  [if  derivative  equations  interlock,  i.e., 

DY(1)  =  function  of  (DY(2)),  DY  (2)  =  function  of  (DY(1)),  these 
equations  must  be  solved,  either  analytically  or  numerically  to 
remove  the  interdependency  before  trying  to  integrate.] 

•  Variables  other  than  T,  Y  and  DY  which  are  used  in  the  expres¬ 
sions  should  be  PAR's  (see  Section  II-B,  STEP  2-2),  or 
variables  defined  in  some  way  within  this  program  segment. 

•  The  exact  number  of  derivatives  to  be  integrated  must  be  indi¬ 
cated  on  the  *RUN  card,  explained  in  Section  VII-A. 

•  To  skip  a  derivative  at  any  time,  simply  set  DY(i)  =  0,  but  make 
certain  that  neq,  specified  on  the  *RUN  card,  corresponds  to  the 
largest  subscript  of  DY  used,  even  though  it  is  desired  to  actually 
integrate  fewer  derivatives. 

The  general  form  of  the  derivative  equation  coding  is 


DY(i)  =  some  function  of  (T,  Y,  DY,PAR, 
constants,  user  variables) 


1 


EXAMPLES: 


0  =  6  +  5.  Ot 


=5»  D  Y  (1 )  =  Y  (1)  +  5.  0*T 


Y  +  bY  =  G^T-2.  0)=S>DY(1)  =  -  B*  Y  (1 )  +  G1  *  (T -2.  0) 


f  =  0  +  2.  Oip 


or  D Y(1 )  =  -PAR(1)*Y(1)  +  PAR(2)*(T-2.  0) 
]  (D  Y(  1 )  =  Y  (2 )  +  2.0*Y(1) 


6  =  0.  5<j>  +  COS(0)  j  |DY(2)  =  0.  5*DY(1)  +  COS(Y(2)) 

Y  +  2{wY  +  w^Y  =  a  +  b  cos  (w^t) 

where  a  =  10.,  b  =  3.,  =  0.05,  £  =  0.  5,  w  =  2. 

^  D  Y  (1 )  =  Y  (2) 

DY(2)  =  10.0  +  3.0*COS(0.05*T)-2.  *0.  5*2.  0*Y(2)-4.  0*Y(1) 


B.  DEFINING  THE  DERIVATIVES  IN  USER  VARIABLES 

Alternatively,  the  user  may  code  his  derivative  equations  using  his  own 
variable  names  if  he  especially  wishes  to  keep  them  more  easily  recognizable 
or  if  he  plans  to  use  vector-matrix  expressions  coded  in  WHELP  variables. 
To  do  so,  however,  he  must  place  within  the  derivative  segment  but  before 
the  derivative  equations,  statements  to  move  each  integrator  output  ( Y (i ) ) 
that  he  plans  to  use  into  his  own  variable  location.  Similarly,  after  his 
equations,  he  must  move  the  values  to  be  integrated  into  the  DY  vector.  This 
may  not  be  done  with  an  EQUIVALENCE  statement,  because  Y  and  DY  are  in 
the  calling  sequence  of  SUBR0UTINE  DERIVS,  and  this  would  be  a  violation 
of  FORTRAN  rules.  It  may  be  done  using  simple  replacement  statements  or 
SUBR0UTINE  M0VE,  whose  calling  sequence  is 

CALL  MOVE  (A,  N,  B) 


where 


A  is  the  first  storage  location  from  which  data  is  to  be  moved. 
It  may  be  specified  as  A  (i),  A  (i,  j),  A(i,  j,k),  or  simply  A 
implying  A(l,  1).  Data  will  be  transferred  by  columns. 

N  is  the  number  of  consecutively  stored  values  tobe  moved. 

B  is  the  first  storage  location  to  which  data  is  to  be  moved, 
specified  in  the  same  way  as  A. 


EXAMPLES: 


i.  e  =  yt  +  6 

y  =  b0(t-5.0) 


*DERIVS 


THETA  =  Y  (1 ) 

THETADT  =  Y(2) 

GAMMA  =  Y  (3) 

THTDTDT  =  GAMMA  *T+THETADT 
GAMMADT  =  PAR(l)*THETA*(T-5.  0) 

DY(1)  =  THETADT 
DY{2)  =  THEDTDT 
DY  (3)  =  GAMMADT 

*ENDDERIVS 


2. 


*DERIVS 


z  =  2 . 0 *z  -  3.  5*^ 

^1  =  ^1T  +  COS 

<p2  =  <pzT  +  sin  ^ 

</> 3  =  +  cos  <fi1 


DIMENSION  PHI(3),  PHIDOT(3) 

C  MOVE  Ys  INT0  USER  VARIABLES. 

Z  =  Y  (1 ) 

CALL  MOVE  (Y(2),3,  PHI) 

C  COMPUTE  DERIVATIVES. 

ZD0T  =  2.0*Z  -  3.  5*PHI(1) 

PHID0T (1)  =  PHI(1)*T  +  C0S  (PHI(2)) 
PHIDOT (2)  =  PHI(2)*T  +  SIN  (PHI(3)) 
PHID0T (3)  =  PHI(3)*T  +  C0S  (PHI(l)) 
C  MOVE  DERIVATIVES  INTO  DYS. 

DY  (1)  =  ZD0T 

CALLM0VE  (PHIDOT,  3,  DY(2)) 
*ENDDERIVS 


3.  A  (3 ,  3 ),  B(3),  C (3 ) 
C  =  A  B 


r> 


*  DERI  VS 

COMMON /USERBLK/ A 
:'TDEC  LARE  A(3,3)  B(3)  CDOT(3) 
CALL  MOVE  ( Y(1 ) ,  3,  B) 

CD0T  =  A*B  $ 

CALL  MOVE  (CD0T,  3,  DY(1)) 


*ENDDERIVS 


$ 


[Note  on  example  3:  Refer  to  Appendix  H  on  WHELP] 
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DEFINING  THE  DERIVATIVES  AS  ENGINEERING 


BLOCKS  (*BLg>CK) 

Since  engineering  systems  are  often  described  oy  block  diagrams  (as  in 
the  example  case,  Section  II- B)  a  special  command  card  *BL$CK  can  be  used 
to  define  the  integration  represented  by  a  block  diagram  so  that  the  user  need 
not  translate  its  contents  into  the  form  Y  =  some  expression. 

A  block  diagram  of  this  form  is  commonly  used  to  represent  a  first 
order  filter  in  an  analog  system: 

6 in  "[jrTTr  -  y(i) 

It  means  simply  that  the  output  of  the  block  equals  the  contents  of  the 
block  multiplied  by  the  input,  or  in  this  case 


y(i>  =  e. 


in  l  s  +  0 


where  1/s  represents  an  integrator.  Solving  this  to  remove  the  s,  we  get 


y(i)  (s  +  0)  =  ein 
y(i)  +  /?*Y(i)  =  e 


y(i)  =  ein  -  £*y(i) 


This  could  now  of  course  simply  be  coded  as 


DY(i)  =  e.n  -0* y(i). 


However,  the  whole  process  of  manipulating  the  variables  (and  thus  possibly 
introducing  errors)  can  be  avoided  by  using  the  special  format: 
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where 

1  denotes  a  first  order  block 

p  is  a  constant  (which  may  be  any  legal  FORTRAN  expression) 

written  without  embedded  blanks 

Y(i)  is  the  dependent  variable  name  for  the  output  of  the  block 

e.  is  a  FORTRAN  expression  for  the  desired  block  input  (It  may 

‘n  have  blanks  within  it,  and  is  terminated  by  the  dollar  sign.  ) 

$  is  the  required  terminator. 

All  items  are  separated  by  blanks. 


Second-order  blocks  can  also  be  specified  directly  by  using  ^BLOCK  2. 
Thus,  a  block  of  this  form 


e .  - 
in 


V  +  *0 

s2  +  PjS  +  P0 

which  is  equivalent  to 

Y.  =  or.e.  +  a  e .  -  p.  Y  .  -  PnY 

i  1  in  0  m  1  l  0  i 


can  be  coded  in  this  format: 


-BLOCK  2  ax  aQ  p:  pQ  Y(i)  Y(j)  e  $ 


where 

2  denotes  a  second-order  block 

q'^.Oq,  Pj,  Pq  are  constants  (which  may  be  any  legal  FORTRAN 
expressions)  written  without  embedded  blanks 

Y(i)  is  the  dependent  variable  name  for  the  output  of  the 

block 

Y(j)  is  another  dependent  variable  name  for  an  inter¬ 

mediate  variable  [it  will  not  be  the  derivative  of  Y ( i )] 

e.  is  a  FORTRAN  expression  for  the  desired  block  input 

(It  may  have  blanks  within  it,  and  is  terminated  by  the 
dollar  sign.  ) 

$  is  the  required  terminator. 

All  items  are  separated  by  blanks. 
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NOTE 


Variables  used  in  *BL0CK  statements  must  follow  the  same 
rules  as  those  used  in  DY(i)  =  ...statements  since  all  *BL0CK 
statements  are  translated  by  the  precompile  program  into  equi¬ 
valent  equations  of  the  form  DY(i)  =  ... 
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SECTION  IV 


INTEGRATION  PACKAGE 

A.  GENERAL  INFORMATION 

The  basic  integration  package  provided  by  ESP  uses  an  algorithm  known 
as  Adams  integration,  which  combines  variable  stepsize  with  variable  order 
integration  to  provide  a  highly  flexible  and  efficient  problem  solving  capa¬ 
bility.  Three  alternative  integration  packages  are  also  available  and  may  be 
selected  at  the  user's  discretion:  they  are  fourth-order  Runge-Kutta,  second- 
order  Runge-Kutta,  and  Hamming  fifth-order  Predictor-Corrector  (which 
uses  fourth-order  Runge-Kutta  as  a  starter).  Each  of  these  alternative 
integrators  may  be  used  with  either  a  fixed  or  variable  stepsize.  (An  expla¬ 
nation  of  the  various  algorithms  used  may  be  found  in  Appendix  E.  ) 

1 .  Options 

Adams  integration  will  be  used  unless  overridden  by  the  user.  To 
select  one  of  the  alternative  packages,  simply  place  the  following  card  before 
all  others  in  your  deck: 


^METHOD  name 

where 

^METHOD  begins  in  column  1 

name  is  RK4  for  fourth-order  Runge-Kutta 

or  RK2  for  second-order  Runge-Kutta 
or  PC  for  Predictor- Corrector 

[See  Appendix  D-3-b  if  you  write  your  own  MAIN  program.] 

The  Adams  integration  package  seems  to  be  highly  successful  with  a 
great  variety  of  problems  and  therefore  should  be  tried  first  for  most  prob¬ 
lems.  Problems  which  require  inputs  at  discrete  intervals,  however,  may  be 
more  suitable  for  Runge-Kutta  integration  (see  Section  IV-C-4). 


4-1 


No  matter  which  integration  algorithm  is  used,  control  of  integration  is 
maintained  by  SUBROUTINE  ESPCTL.  It  computes  the  initial  stepsize  (if 
variable),  calls  the  specified  integration  routine,  calls  the  switching  routines, 
checks  for  occurrence  of  switches,  and  handles  printing  and  storage  of  plot 
data. 

2.  Stepsize  Selection 

Immediately  after  the  first  call  to  SUBROUTINE  DERIVS,  ESPCTL  will 
determine  the  initial  H  in  one  of  the  following  ways: 

•  If  both  Y.  and  Y.  are  nonzero,  then 

1  1  ’ 

H  =  min  Y./Y.l  *0.2^ 

•  If  either  Y.  or  Y.  is  zero  for  all  i,  then 

l  l  ’ 

H  '  5l2~ 

•  The  user  may  define  H  in  ICCOMP,  and  this  H  will  be  tried 
first  (see  below). 

•  The  user  may  select  fixed  stepsize  integration  (with  RK2,  RK4, 
or  PC  only)  by  setting  FEXSTP  >  0  (see  below). 

After  initial  stepsize  selection,  stepsize  control  proceeds  differently  in 
each  integration  package  and  is  documented  in  the  description  of  each  method 
given  below. 

3.  User  Control  of  Stepsize 

In  order  to  change  any  of  the  stepsize  variables,  the  user  should  place 
the  following  common  block  in  ICCOMP  and  use  arithmetic  statements  to 
define  those  variables  he  wishes  to  change: 
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■>  < 


where 


HP 


H 


F1XSTP 

HMIN 


HMAX 


is  the  current  printing  interval.  This  is  normally  changed 
from  the  *RUN  card  but  may  be  changed  by  the  user's 
program  during  the  run  and  must  be  >  0. 

will  be  the  initial  stepsize  tried,  and  the  current  stepsize 
during  execution.  It  may  be  set  only  in  ICCOMP  or  at 
switching  times  and  must  satisfy  HMIN  ^  H  <  HMAX. 

(default  =  0)  if  set  >  0,  causes  H  =  FIXSTP  at  all  times. 

(default  =  0)  the  lower  limit  on  the  stepsize.  HMIN  >  0 
causes  printing  of  a  warning  message  and  continuation  of 
integration  with  acceptance  of  errors  whenever  H  <  2.0  X 
HMIN. 

(default  =  1.  0E50)  the  maximum  stepsize  permitted. 


4.  Controlling  Solution  Accuracy 

Solution  accuracy  may  be  controlled  by  the  user  through  either  or  both 
of  two  variables,  which  may  be  input  on  run-time  data  cards.  The  exact  use 
of  these  variables  depends  upon  the  integration  algorithm  selected  and  is 
explained  in  the  discussion  of  each  given  below.  The  cards  should  be  placed 
after  all  derivative,  input,  and  output  coding  but  before  the  *RUN  card.  The 
formats  are: 


*EPS 

€ 

*Q 

*1 

q2  *  *  *  qn 

$ 

where 


Default  value  of  €  =  1.  E  -  6 
Default  value  of  =  l.E  -  10 
*EPS  and  *Q  start  in  column  1 

q.  can  be  either  =  a  constant  alone  (assumed  to  be  control  for 

next  variable) 

or  Qi  =  constant 

or  ALL  =  constant 

$  is  a  required  terminator  for  *Q 
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EXAMPLES: 

*EPS  l.E-2 

*Q  Q2  =  l.E-20 

l.E-8 

Q5  =  1.E5 

$ 

*Q  ALL  =  l.E-6 

$ 

B.  ADAMS  INTEGRATION 

The  Adams  integration  package  features  both  variable  stepsize  and 
variable  order  integration.  Since  it  is  able  to  dynamically  adjust  both  the 
stepsize  and  the  order  according  to  the  success  of  each  integration  step,  it  is 
capable  of  solving  a  wide  diversity  of  problems  with  both  accuracy  and  speed. 

On  the  first  call  to  ADAMS,  an  appropriate  stepsize  H  is  computed 
using  the  H  supplied  by  ESPCTL  as  a  starting  point,  a  flag  IFAIL  is  set  to  0, 
the  order  K  is  set  to  1,  and  all  necessary  coefficients  of  formulas  are 
initialized  and  computed.  On  subsequent  calls,  IFAIL  is  reset  to  0  and  H 
retains  the  value  assigned  in  the  previous  call  and  is  merely  tested  to  ensure 
that  it  is  within  the  precision  limits  of  the  computer.  Coefficients  needed  for 
integration  are  then  recomputed  only  if  H  has  changed. 

From  this  point  on,  the  sequence  of  events  is  the  same  for  each  call  to 
ADAMS.  A  solution  is  predicted  and  DERIVS  is  called  to  evaluate  the  deriva¬ 
tives  at  the  predicted  solution.  The  local  error  is  estimated  at  order  K,  K-l, 
and  K-2,  and  if  necessary  the  order  is  lowered  for  the  next  step. 

If  the  errors  are  within  an  acceptable  range,  the  step  is  considered 
successful,  the  predicted  solution  is  corrected  and  the  derivatives  re-evaluated 
at  T  +  H.  Differences  are  updated  for  the  next  step,  and  the  best  order  and 
stepsize  are  determined  for  the  next  step  before  control  is  returned  to  ESPCTL. 

If  the  errors  are  not  acceptable,  T  is  reset  to  T-H,  the  flag  IFAIL  is 
incremented  by  1  and  then  tested.  On  the  first  and  second  failures,  H  is 
halved  and  the  order  retained  before  retrying  integration;  on  the  third  failure. 
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the  order  K  is  also  reduced  to  1.  If  more  than  three  failures  occur,  an 
optimum  H  is  also  computed  before  retrying.  Again  the  size  of  H  is  tested, 
and  if  too' small  for  machine  precision,  the  error  tolerance  is  doubled,  and 
KFLAG  is  set  to  0  so  that  a  warning  message  will  be  printed: 


"REQUESTED  ACCURACY  NOT  ACHIEVED  AT  T  =. 
REMAINDER  OF  SOLUTION  IS  SUSPECT.  " 


The  error  control  method  of  ADAMS  utilizes  both  stepsize  and  order 
variation  to  keep 


r  ,  1/2 

neq 

* 


where 


E^  -  estimate  of  the  error  in  Y^  mode  in  the  current  step 

Q.  *  maximum  thus  far  of  Q.  and  |Y.  |  (updated  at  the  end  of  each 
integration  step) 


(See  Ref.  6. ) 


Bmn 


^Bi 


C.  RUNGE-KUTTA  INTEGRATION 

Both  fourth-order  Runge-Kutta  and  second-order  Runge  may  be  used 
with  either  a  fixed  stepsize  or  a  variable  stepsize.  Each  is  implemented  in 
its  own  subroutine,  but  these  routines  are  parallel  in  logic  and  sequence,  and 
differ  only  in  the  integration  equations  and  constants  used.  The  following 
paragraphs  refer  to  fourth-order  Runge-Kutta  with  differences  applying  to 
second-order  Runge-Kutta  noted  in  parentheses. 

1 .  Fixed  Stepsize 

If  the  variable  FIXSTP  is  set  >  0  by  the  user  (see  Section  IV-A-3),  the 
stepsize  H  =  FIXSTP  at  all  times,  and  no  error  testing  of  any  kind  is  done. 
Each  call  to  the  integration  routine  causes  two  single  integration  steps,  namely 
from  T  to  T+H  and  from  T+H  to  T+2H. 

2.  Variable  Stepsize 

If  FIXSTP  $  0  (default  is  0. ),  integration  begins  with  the  stepsize  H  as 
determined  in  ESPCTL.  A  double  step  is  taken,  from  T  to  T+2H,  and  then 
two  single  steps  are  taken,  from  T  to  T+H  and  from  T+H  to  T+2H.  The  total 
error  is  computed  and  compared  with  the  permitted  error  bounds.  If  not 
acceptable,  stepsize  doubling  is  prevented  and  H  is  compared  with  HMIN  to 
determine  if  the  stepsize  can  be  halved.  If  not,  integration  continues  using 
the  current  H  and  a  warning  message  is  printed: 

"REQUESTED  ACCURACY  NOT  ACHIEVED  AT  T  = _ 

REMAINDER  OF  SOLUTION  IS  SUSPECT." 

If  H  can  be  halved,  it  is  and  the  above  process  is  repeated. 

If  the  errors  are  within  acceptable  limits,  they  are  further  tested.  If 
less  than  0.  5%  (1%  for  RK2)  of  the  error  bounds,  the  stepsize  H  is  permitted 
to  double  for  the  next  integration  step. 
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3. 


Error  Control 


The  error  aLlowed  in  the  computation  of  Y (i ) s  is  controlled  by  requiring 

that 

1=1 

and  Q(I)  is  initially  set  to  MAX(Q(I),  |  Y0(I)|)  and  then  continuously  updated  to 

Q(I)  =  MAX(Q(I),  | Y(I)|  ,  |  Y (I)  |) 

(T  =  T-H) 


The  default  value  of  EPS  is  l.E-6  and  of  Q(I)  is  l.E-10,  but  the  user 
may  change  them  to  suit  his  problem.  (See  Section  IV-A-4.  ) 

4.  Inputting  Values  at  Discrete  Intervals 

Since  the  user  sometimes  wishes  to  input  noise  or  compute  values  at 
discrete  and  predictable  intervals  during  integration  and  because  the  number 
of  evaluations  of  the  derivatives  is  different  for  each  integration  routine,  a 
special  flag,  FIRSTP,  which  signals  the  beginning  of  each  integration  step 
(or  pair  of  steps)  has  been  added.  To  use  this  flag,  include  the  following  card 
in  the  derivative  segment  of  couiug 

C0MM0N/RKC0NT/ FIRSTP 


and  ,est  FIRSTP  to  determine  when  to  input  values.  FIRSTP  =  0.  normally, 
but  is  set  to  1.  0  by  the  integration  routine  at  the  beginning  of  each  step  in  the 
fixed  step  mode  or  at  the  beginning  of  each  pair  of  steps  in  the  variable  step 
mode. 
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Simulations  with  noise  may  be  easily  set  up  by  using  ^METHOD  RK4  or 
RK2,  setting  FIXSTP  *  0.  (See  Section  IV-A-3)  and  inputting  a  noise  value 
whenever  FIRSTP  =1.0  (See  the  example  below.). 

Alternatively,  the  automatic  stepsize  control  feature  can  be  retained 
and  a  crude  kind  of  switching  capability  achieved  by  using  ^METHOD  RK4 
or  RK2  with  a  variable  stepsize  and  simply  testing  the  FIRSTP  flag  (as  shown 
in  the  example  below).  This  will  permit  successful  introduction  of  discrete 
or  discontinuous  values  at  the  beginning  of  each  integration  step,  irrespective 
of  the  step  size. 

EXAMPLE: 

*  METHOD  RK4 
*DERIVS 

C0MM0N/RKC0NT/ FIRSTP 

IF  (FIRSTP  .  NE  1)  GO  TO  5 

AN0ISE  =  RANF(O) 

Y  (1)  =  Y(l)  +  AN0ISE 
5  CONTINUE 

D  Y  (1 )  =  .  .  . 

*ENDDERIVS 
*ICC0MP 

COMMON/ST PCON/HP,  H,  FIXSTP,  HMIN,  HMAX 

FIXSTP  =  0.  02 


*ENDIC 


D.  PREDICTOR-CORRECTOR  INTEGRATION 

Hamming's  fifth-order  Predictor-Corrector  is  the  algorithm  used,  but 
since  it  is  not  se if- starting ,  fourth-order  Runge-Kutta  is  used  to  start  the 
solution  at  T^  and  to  restart  the  solution  after  discontinuities  or  difficulties 
are  encountered. 

Using  the  stepsize  H  as  determined  in  ESPCTL,  steps  1-3  are  taken 
with  fourth-order  Runge-Kutta  and  the  errors  are  checked  at  T^.  If  the  error 
on  this  step  exceeds  the  error  bounds,  the  stepsize  is  halved  and  the  solution 
is  restarted  from  Tq.  If  the  error  is  acceptable,  step  4  is  taken  with  Runge- 
Kutta  and  step  5  with  Predictor  -  Corre  ctor. 


Fig,  1.  Step  Sequence  for  Starting  Procedure 


1 .  Variable  Stepsize 

Once  the  solution  has  been  started  in  this  way,  error  checks  are  made 
continuously,  and  the  stepsize  is  halved  when  the  error  exceeds  the  bounds 
and  permitted  to  double  for  the  netft  step  when  the  error  is  less  than  1%  of 
the  bounds.  To  prevent  excessive  interval  halving  if  the  problem  happens  to 
be  ill-conditioned,  a  counter  is  used  to  monitor  stepsize  halvings.  Each  time 
the  solution  is  restarted  with  Runge-Kutta,  the  counter  (KC<ftUNT)  is  started 
at  zero.  Each  time  the  stepsize  is  decreased,  KC$UNT  is  decreased  by  one, 
and  after  each  successful  step  it  is  increased  by  one,  but  never  permitted  to 
exceed  zero.  If  KCOUNT  becomes  less  than  -4,  a  warning  message  is  printed: 
"SOLUTION  APPEARS  ILL-CONDITIONED  AND  IS  BEING  RESTARTED.  THE 
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FOLLOWING  Y'S  EXCEED  THE  ERROR  BOUND",  followed  by  the  Y(I)s.  T  is 
reset  to  T-H,  and  the  solution  is  restarted  using  Runge-Kutta. 


The  above  description  of  stepsize  and  error  control  assumes  that  the 
variable  HMIN  has  its  default  value  of  zero.  If,  however,  HMIN  is  set 
greater  than  zero  (See  Section  IV-A-3)  and  the  error  is  too  large,  the  step- 
size  H  is  halved  if  possible  and  the  solution  continued.  But,  if  it  cannot  be 
halved  without  making  it  less  than  HMIN,  it  retains  its  value,  the  solution 
continues  and  a  warning  message  is  printed:  "REQUESTED  ACCURACY  NOT 

ACHIEVED  AT  T  = _ ,  REMAINDER  OF  SOLUTION  IS  SUSPECT.  "  This 

option,  it  may  be  seen,  may  produce  less  accuracy  but  in  some  cases  more 
speed.  The  user  is  advised  to  consider  any  warning  messages  he  receives 
and  to  base  his  selection  of  HMIN  on  the  nature  of  his  problem  and  the  desired 
results. 


2.  Fixed  Stepsize 

Although  the  Predictor-Corrector  integration  package  is  intended  for 
use  as  a  variable  stepsize  method,  it  can  also  be  used  with  a  fixed  stepsize  by 
setting  FIXSTP  >  0.  (See  Section  IV-A-3.  ).  In  this  case,  all  error  checking 
is  skipped  and  no  interval  halving  or  doubling  occurs. 

3.  Error  Control 


The  error  allowed  in  the  computation  of  Y(i)s  is  controlled  by  requiring 

that 


EPS  > 


|  ERR(/)R(I)| 
Q(I) 


for  all  I 


and  Q(I)  is  initially  set  to  MAX(Q(I),  Y 0 ( I)  | )  and  then  continuously  updated  to 


Q(I)  =  MAX(Q(I),  |  Y(I)|) 


The  default  value  of  EPS  is  l.E-6  and  of  Q(I)  is  l.E-10,  but  the  user 
may  change  them  to  suit  his  problem  (See  Section  IV- A -4). 
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SECTION  V 


DISCONTINUITIES 


In  programming  a  system  of  differential  equations,  the  need  frequently 
arises  for  a  means  to  model  accurately  various  types  of  discontinuities  and 
nonlinearities  which  are  part  of  the  system.  Therefore,  several  special 
features  have  been  built  into  ESP  to  handle  the  most  common  types  of  non- 
linearities,  and  with  the  aid  of  a  little  imagination  (some  examples  will  be 
given)  nearly  any  desired  characteristic  can  be  produced  by  modifying  one  of 
the  features. 

It  is  important  to  understand  that  these  features  are  provided  not  merely 
for  convenience,  however!  Since  the  integration  algorithms  available  with 
ESP  by  and  large  assume  that  they  are  working  on  continuous  and  reasonably 
well-behaved  derivatives,  haphazard  introduction  of  discontinuities  by  the 
user  can  cause  enormous  problems  and  errors.  The  user  is  strongly  urged 
to  be  certain  that  discontinuities  are  introduced  only  by  means  of  one  of  the 
devices  documented  below  and  that  the  constraints  mentioned  with  regard  to 
their  use  be  closely  observed.  (Section  IV-C-4  discusses  one  other  method 
of  introducing  discontinuities,  which  may  be  used  with  Runge-Kutta  integra¬ 
tion.  ) 

These  special  features,  which  may  be  used  only  in  the  *DERIVS  seg¬ 
ment,  consist  of  SWTCH's,  which  detect  sign  changes  in  an  expression  and 
restart  integration,  SWMEM's,  which  represent  hysteresis  nonlinearities 
and  also  restart  integration,  and  an  EVENT  locator,  which  detects  and  reports 
the  occurrence  and  timing  of  any  user-specified  event  but  which  does  not 
affect  the  integration  process.  Basic  use  of  these  features,  which  is  fairly 
straightforward,  will  be  documented  first,  and  the  later  part  of  this  section 
will  be  devoted  to  extended  usage,  a  description  of  how  the  switches  work, 
and  some  details  and  considerations  regarding  timing  and  accuracy. 


♦ 
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A. 


DETECTING  A  SIGN  CHANGE  (*SWTCH) 


The  *SWTCH  command  detects  a  change  of  sign  in  its  control  or  input 
expression,  locates  the  time  of  this  change  within  a  specified  degree  of  accu¬ 
racy  (see  Section  V-F),  assigns  itself  an  output  according  to  the  sign  of  the 
input  expression,  produces  an  automatic  print  point  at  the  switch  time  (which 
the  user  may  suppress),  and  restarts  the  integration  from  the  switch  time. 

It  is  useful,  therefore,  in  producing  an  accurate  discontinuous  driving  func¬ 
tion  to  a  derivative  equation  and  in  permitting  the  user  to  detect  the  exact 
time  of  a  switch  and,  if  he  wishes,  to  perform  some  specific  act  at  that  time. 
Other  possible  uses  will  be  illustrated  in  the  examples. 

The  general  form  of  the  *SWTCH  command  is 


*SWTCH  i  0+  $  0  $  control  expression^  $ 


where 

*SWTCH  starts  in  column  1 

i  (1  <  i  <  50)  is  the  number  of  this  switch 

0  +  is  any  legal  F0RTRAN  expression  which  will  be  the  output  if 
the  control  expression  >  0. 

0  is  any  legal  F0RTRAN  expression,  which  will  be  the  output  if 
the  control  expression  ^  0. 

control  expression.  is  any  legal  FORTRAN  expression  involving 

only  T,  Y,  PAR,  system  functions  and  constants. 
(For  use  of  other  variables,  see  Section  V-D-l.  ) 

$  is  a  required  terminator  of  the  0+,  0  ,  and  control  expressions 

The  ESP  precompiler  breaks  up  the  *SWTCH  card  coding  into  the  input 
(control  expression),  which  it  writes  as  part  of  SUBR0UTINE  SWINPT  in  the 
form  VALUES(i)  =  control  expression^,  and  the  output  computation  which  it 
writes  as  part  of  SUBR0UTINE  DERIVS. 


There  are  two  output  variables  available  to  the  user  resulting  from 
the  *SWTCH  statement.  The  first,  SWCHi,  is  available  only  in  the  deriva¬ 
tive  segment  of  coding.  To  use  it  elsewhere,  such  as  in  0UTPUT,  the  user 
must  compute  it  himself  (see  Section  V-D-2).  The  second,  SWTCH(i),  is 
available  in  DERIVS,  0UTPUT,  SWINPT,  SWMEMN  and  in  any  other  routine 
in  which  the  common  block  SWTCHS  appears.  The  variables  contain  the  fol¬ 
lowing  information: 


SWCHi 


0+  if  control  expression^  >  0. 
0  if  control  expression^  <  0. 


SWTCH(i)  where:  |SWTCH(i)|  is  one  larger  than  the  number  of  sign 

changes  made  thus  far  by  the  control  expression, 
and  is  normally  a  floating  point  integer 

The  sign  of  SWTCH(i)  is  the  current  sign  of  the 
control  expression 

SWTCH(i)  serves  as  a  signal  to  the  user  that  a 
switch  has  just  occurred:  On  the  first  call  to 
DERIVS  following  a  switching,  each  SWTCH(i) 
which  has  been  toggled  has  its  absolute  value 
increased  by  1.  5.  See  example  2  below  for  a 
way  to  utilize  this  trait,  and  see  Section  V-E  for 
detail  on  the  exact  sequence  of  events  when  a 
switch  occurs. 


EXAMPLES: 

1.  The  example  problem  in  Section  II  shows  a  typical  use  of  *SWTCH: 


Input - 


o 

• 

L 

1 

Output 


This  is  coded  as 


*SWTCH  1 


1.0  $ 


0.  0  $  Input  $ 
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which  produces  the  following  results: 


If  input  <  0.,  SWCH1  =  0.,  SWTCH(l)  =  -  N  (N  =  number  of  sign 

changes  +  1) 

If  input  >  0.,  SWCH1  =  1.  0,  SWTCH(l)  =  +  N 

Also,  when  the  switch  is  detected  and  located  within  HSW(l)  of  the 
time  it  occurs,  then 

SWCH1  =  0.,  if  input  <  0. 

=  1 .,  if  input  >  0. 

and  SWTCH(l)  =  SIGN ( | N  |  +  1 .  5,  Input) 

In  this  example,  SWCH1  is  the  relevant  output  and  is  used  as  a  term 
in  th'’  expression  for  DY(1) 

DY(’)  =  G1+SWCH1  -  B*Y(1) 

2.  To  detect  the  exact  time  at  which  subroutines  are  to  be  called  to  re¬ 
define  a  number  of  program  constants,  the  following  arrangement 
could  be  used.  Notice  that  no  value  is  assigned  toSWCHl  because 
the  output  variable  of  interest  here  is  SWTCH(i),  and  that  the  coding 
makes  use  of  the  fact  that  SWTCH(i)’s  are  nonintegers  exactly  at 
switch  times. 

(assume  PAR(l)  =time^,  PAR(2)  =  time^,  etc.) 

*SWTCH  1  0.  $  0.  $  T-PAR{1)  $ 

IF(SWTCH(1)  .NE.  AINT (SWTCH(1)))CA LL  DUMDUM  1 
*SWTCH  2  0.  $  0.  $  T-PAR(2)  $ 

IF(SWTCH(2)  .NE.  AINT(SWTCH(2)))CALL  DUMDUM2 

3.  To  produce  a  sample  and  hold  at  times  t^,  t^,  t^,  ...  a  similar  but  more 
abbreviated  setup  can  be  used,  employing  only  one  switch  which  will 
detect  a  sign  change  as  each  successive  time  is  reached.  Dimension 

a  vector  TSAMP  of  length  N  and  store  the  desired  sample  times  T., 
^2’  *  *  ’  ^ m  Into  it.  Initialize  SAMPLE,  PAR(3)  =  t^  and  1=1  and 
then  use  the  following  coding: 

+SWTCH  4  0.  $  0.  $  T  -PAR  (3)  $ 

IF(SWTCH{4)  .EQ.  AINT  (SWTCH(4)))  G0  T0  5 
C  (THIS  SECTION  WILL  BE  EXECUTED  ONLY  AT  SWITCH  TIMES.  ) 
SAMPLE  =  .  . . 

1  =  1+1 

PAR  (3)  =  TSAMP(I) 

5  CONTINUE 


c 
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In  this  example  it  is  assumed  that  t_  <  t,  <  t_  .  . .  <  t  .  As  execution 

0  l  2  n 

begins,  this  switch  is  first  toggled  when  t  equals  t^  and  at  that  time 
SWTCH(4)  has  1.  5  added  to  it  so  that  the  IF  test  will  fail.  Thus,  SAMPLE 
is  computed  and  PAR(3)  is  reset  to  the  next  sample  time.  Since  immediately 
after  this  the  .  5  is  stripped  from  SWTCH(4)  this  coding  will  be  then  bypassed 
until  t  reaches  the  new  value  of  PAR(3).  Note  that  this  represents  one  of  the 
few  cases  in  which  it  is  permissible  to  store  a  time -dependent  value  in  PAR 
and  to  change  the  input  to  a  switch  in  a  discontinuous  manner. 


B.  HYSTERESIS  NON  LINEARITIES  (*SWMEM) 

A  hysteresis  nonlinearity,  of  the  general  type  illustrated  in  Fig.  2, 
can  be  modeled  using  the  *SWMEM  cards  explained  below. 


INPUT 


It  will  determine  the  proper  location  and  output  of  the  hysteresis  within  the 
required  accuracy  (see  Section  V-F),  inform  the  user  by  setting  a  flag  when¬ 
ever  a  discontinuity  occurs,  indicate  whether  the  output  is  in  the  linear,  dead¬ 
band  or  saturation  regions,  and  restart  the  integration  at  each  discontinuity. 
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There  are  three  special  control  cards  which  may  be  used  to  implement 
this  option: 


*SWMEM 


(required)  defines  input  to  a  SWMEM. 


;:SWMEMDATA  (optional)  defines  the  characteristics  of  the  SWMEM. 


*SWMEMSET 


(optional)  initializes  in  saturation  instead  of  at  zero. 


1.  Defining  Input  to  a  SWMEM 

The  input  to  the  hysteresis  is  defined  on  the  *SWMEM  card  which  is 
placed  in  the  derivative  segment  of  coding.  The  ESP  precompiler  will  write 
the  input  as  part  of  SUBROUTINE  SWMEMN  in  the  form  VALUES(i)  =  input.. 
The  format  is 


*SWMEM 


i  inpub  $ 


where 


*SWMEM  starts  in  column  1 

i  (1  £  i  £  50)  is  the  number  of  this  SWMEM 

inpub  is  any  legal  FORTRAN  expression  involving  only  the  variables 
T,  Y,  PAR,  system  functions  and  constants.  (See  Section 
V-D-l  if  other  variables  must  be  used.  ) 

$  is  a  required  terminator 


Z3 


2.  Defining  Output  of  a  SWMEM 

There  are  two  separate  output  variables  from  SWMEM's,  parallel  in 
nature  to  those  from  SWTCH's.  The  first,  SWMi,  is  automatically  com¬ 
puted  and  made  available  to  the  user  in  the  derivative  segment.  To  use  it 
elsewhere,  the  user  must  compute  it  himself  (refer  to  Section  V-D-2).  The 
second  output  variable,  SWMEMfi,  4),  is  available  in  any  routine  where  the 
common  block  SWTCHES  appears.  The  variables  contain  the  following 
information: 

th 

SWMi  is  the  actual  output  value  of  the  i  SWMEM 

SWMEM(i,  4)  is  normally  a  floating  point  integer  indicating  the 
present  position  on  the  hysteresis  by  its  value: 

-2.0  indicates  negative  saturation 
-1.0  indicates  slope  on  negative  side  ' 

0.  0  indicates  deadband  ^ 

1.  0  indicates  slope  on  positive  side 

2.  0  indicates  positive  saturation 

SWMEM (i,  4)  also  signals  the  user  that  a  "corner  has 
just  been  turned"  on  the  hysteresis:  On  the  first  call 
to  DERIVS  following  a  SWMEM  discontinuity,  the 
absolute  value  of  each  SWMEM(i,  4)  which  has  changed 
state  is  increased  by  0.  5.  After  the  derivative  equa¬ 
tions  are  evaluated,  the  SWTCH's  and  SWMEM's  are 
reevaluated  and  the  0.  5  removed  before  the  integration 
is  restarted.  (This  signal  may  be  tested  and  used  in  the 
same  ways  that  SWTCH(i)'s  are  used  in  the  examples 
(refer  to  Section  V-E  for  more  detail  on  the  exact 
sequence  of  events.  ) 

3.  Defining  the  Characteristics  of  a  SWMEM 

Generally  the  user  will  want  to  define  the  constants  Cl  through  CIO 
(see  Fig.  3)  characterizing  his  SWMEM,  although  they  do  have  default  values 
for  the  simplest  case.  Constants  C3,  C4,  C8,  and  C9  are  the  slopes.  How¬ 
ever,  an  infinite  slope  is  defined  by  setting  the  corresponding  C^  equal  to  zero. 


If  the  user  wishes  to  know  what  path  he  is  following  on  the  hysteresis,  he 
may  store  the  past  value  of  SWMEM(i,  4)  so  he  will  know  where  he  is 
coming  from  at  each  "corner." 
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^1 


Fig.  3.  SWMEM  Characteristics 

The  C's  are  usually  defined  as  part  of  the  run-time  data  cards, 
meaning  that  they  are  placed  somewhere  between  *ENDIC  (if  used)  and 
*RUN,  and  are  picked  up  in  the  same  way  that  PAR's  are  picked  up  from 
*PAR  cards.  The  format  for  inputting  C's  is 


*SWMEMDATA 


10 


'10 


$ 

$ 


whe  re 


*SWMEMDATA  starts  in  column  1 

i  is  the  number  of  the  corresponding  SWMEM 

c.  is  either  a  constant  alone  or  Cj  =  constant,  and  blanks  are 

■'  separators.  If  a  constant  appears  alone,  it  is  assumed  to 

be  the  value  of  the  next  c:.  If  no  value  is  given  for  c-,  its 
default  value  will  be  usea.  " 

$  is  a  required  terminator 


I 
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In  general  Cl  must  be  >  C6,  except  for  special  case  1  below. 


Default  values  are  C5  =  1,  CIO  =  -1,  all  other  Cs  =  0.,  giving 


Out 

_ 1 

put 

O 

• 

I 

SPECIAL  CASES: 

1.  Cl  =  C2  =  C6  =  C7  is  permitted,  giving  something  like  this: 


2.  For  a  symmetric  nonlinearity  (corresponding  to  forming  quadrant  III 
by  rotating  quadrant  I  through  180  degrees  about  the  origin),  only  Cl 
through  C5  need  be  defined  with  the  result  that 

C6  =  -Cl  C9  =  C4 

C7  =  -C2  C8  =  C3 

CIO  =  -C5 

3.  If  it  is  necessary  to  compute  any  of  the  C  values,  they  may  be  defined 
as  C0NSTS  in  ICC0MP,  if  the  following  common  block  is  added: 


C0MM0N/SWHPAR/NCHNG,  NALTER,  ISWTYP,  KSV, 
C0NSTS(5O,  10),  SWSET (50) 

C  values  defined  in  ICC0MP  will  supersede  any  that  are  input  on 
#SWMEMDATA  cards,  but  no  error  check  will  be  made  on  them. 
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4.  Initializing  a  SWMEM 

Normally,  the  function  is  initialized  on  the  region  corresponding  to 
starting  from  zero,  but  it  may  be  initialized  on  the  region  corresponding  to 
starting  from  saturation  by  using  the  *SWMEMSET  card  among  the  run-time 
cards.  The  format  is 


where 


*SWMEMSET  nL  n2 


*SWMEMSET  starts  in  column  1 


m  is  the  number  of  the  *SWMEM  to  be  set  and  blanks  are  separators 

$  is  a  required  terminator  of  the  list 

The  m  are  not  retained  from  run  to  run,  so  it  is  necessary  to  redefine  them 
for  different  runs. 


C.  LOCATING  EVENTS  WHICH  DO  NOT  AFFECT  INTEGRATION 

The  EVENT  capability  is  useful  for  finding  events  which  do  not  intro¬ 
duce  discontinuities  into  the  differential  equations.  It  detects  the  time  of 
occurrence  of  any  event  specified  in  the  user-written  SUBROUTINE  EVENTS 
within  the  timing  accuracy  specified  on  the  *HSWE  card  (see  Section  V-F). 
The  event  is  recognized  as  a  change  of  sign  in  the  input  or  event -defining 
expression,  and  any  event  can  therefore  occur  many  times  within  a  run. 
Having  located  an  event,  ESP  then  interpolates  T  and  Y  to  the  event  time  and 
calls  the  user- supplied  SUBROUTINE  NOTIFY,  as  its  only  response  to  the 
event.  It  does  no  printout  and  it  does  not  in  any  way  affect  the  integration 
process  or  the  rest  of  the  run. 
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To  use  the  EVENT  locator,  three  things  must  be  done: 


•  Place  an  *NEVENT  card  among  the  run-time  data  cards: 


*NE VENT  n 


where 

*NEVENT  starts  in  column  1 
n  is  the  number  of  events  to  be  defined 

•  Supply  the  subroutine  to  define  the  functions  which  determine  the 
e  ve  nt  s : 

SUBR0UTINE  E VENTS(VA  LUES,  T,  Y) 

DIMENSION  VALUES  (1),  Y(l) 

[Other  common  and  dimension  statements  as  needed] 

VALUES(l)  =  expression  determining  event^ 

• 

VALUES(n)  a  expression  determining  event 
RETURN  n 

END 

•  Supply  the  subroutine  to  receive  notification  that  event  number  IEVENT 
has  occurred  at  the  given  values  of  T  and  Y.  This  routine  will  be 
called  once  for  each  event  occurrence,  in  the  order  in  which  events 
occur.  Its  format  is; 

SUBROUTINE  NOTIFY (T,  Y,  IEVENT) 

DIMENSION  Y(1 ) 

[Other  common  and  dimension  statements  as  needed] 

[Statements  defining  how  EVENT  information  is  to  be  used] 

RETURN 

END 


EXAMPLE: 

SUBROUTINE  NOTIFY(T,  Y,  IEVENT) 

DIMENSION  Y(1 ) 

PRINT  100,  IEVENT,  T 

100  F0RMAT(1HO,  *EVENT*,  14,  ^OCCURRED  AT*,  E8.  2) 
RETURN 
END 
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D.  SWTCH's  and  SWMEM's;  EXTENDED  USAGE 

Frequently  the  standard  usage  of  SWTCH's  and  SWMEM's  is  too  con¬ 
fining  for  the  user's  needs,  either  because  of  the  limitations  on  how  inputs 
may  be  defined  or  because  the  desired  output  values  are  not  automatically 
available  in  all  routines.  The  following  paragraphs  illustrate  several  ways 
that  these  limitations  on  both  input  and  output  can  be  circumvented. 

1 .  Alternate  Ways  to  Define  SWTCH  and  SWMEM  Inputs 

Basic  usage  of  *SWTCH  and  *SWMEM  limits  the  form  of  their  inputs  to 
simple  FORTRAN  expressions  using  T,  Y,  PAR,  system  functions  and 
constants  only,  because  of  the  way  these  statements  are  translated  by  the 
precompiler  into  the  input  routines  SWINPT  and  SWMEMN.  Since  it  is 
sometimes  necessary  either  to  use  other  variables  or  to  execute  a  series  of 
statements  to  define  a  switch  input,  alternate  means  of  defining  switch  inputs 
are  available.  The  first,  and  probably  simplest,  is  for  the  user  to  write 
a  function  subprogram  which  defines  his  input;  the  second  is  to  simply  write 
the  entire  SWINPT  or  SWMEMN  routine  himself.  In  either  case  the  user 
should  remember  that  T  and  Y  are  always  updated  for  the  purpose  of  com¬ 
puting  switch  inputs,  and  that  any  time -dependent  variables  used  to  compute 
switch  inputs  should  themselves  be  computed  within  the  function  subprogram 
or  within  the  user-written  SWINPT  (SWMEMN)  so  they  too  will  be  properly 
updated. 

a.  User-Written  Functions 

The  major  advantage  in  using  function  subprograms  to  define  switch 
inputs  is  that  it  permits  the  user  to  combine  the  convenience  of  the 
*SWTCH(*SWMEM)  card  with  almost  total  flexibility  in  defining  his  input. 

He  may  use  common  blocks  to  pass  his  own  variables  to  the  switch  input 
computation,  and  he  may  use  as  many  F0RTRAN  statements  as  he  wishes  to 
define  his  actual  input.  The  following  example  illustrates  some  possibilities 
of  this  approach: 


5-12 


EXAMPLE: 


If  the  switch  input  is 

if  K1  sr  1,  input  =  Y (3 )  *  10.  -BETA 
if  K1  =  0,  input  =  Y{3)  *B ETA  +  Y(2) 

and  the  desired  output  is 

if  input  >  0,  output  =  1.0 
if  input  £  0,  output  =  0.0 

the  switch  can  be  coded  by  putting  the  following  cards  in  the  derivative 
section 

C0MM0N/BL0CK1/K1,  BETA 

• 

*SWTCH  1  1.0  $  0.0  $  SWFUNC(Y)  $ 

and  writing  the  following  function  subprogram  to  be  placed  before  the 
derivative  coding 

FUNCTION  SWFUNC(Y) 

C0MM0N/BL0CK1/K1,  BETA 
DIMENSION  Y(1 ) 

IF  (K1  .  EQ.  0)  G0  T0  5 
SWFUNC  =  Y(3)  *  10.  -BETA 
RETURN 
5  CONTINUE 

SWFUNC  =  Y(3 )  *BETA  +  Y(2) 

RETURN 

END 


b. 


User-Written  SWINPT  and  SWMEMN 


If  the  user  has  many  switch  inputs  which  require  extensive  computa¬ 
tion,  he  may  prefer  to  simply  write  his  own  input  routines  rather  than  write 
many  functions.  As  with  the  function  subprograms,  any  number  of  common 
blocks  and  computations  may  be  included  when  the  user  writes  his  own  sub¬ 
routines  SWINPT  or  SWMEMN.  He  must,  however,  write  these  routines  in 
the  form  in  which  ESP  expects  them,  as  explained  below,  and  be  careful  to 
place  them  after  his  job  control  cards  but  before  the  *DERIVS  segment, 
which  will  result  in  their  being  used  instead  of  the  dummy  routines  written  by 
the  precompiler.  Either  routine  or  both  may  be  user-written,  but  if  the  user 
writes  his  own  routine  for  SWTCH  inputs  (SWMEM  inputs),  he  must  define  all 
of  his  SWTCH  inputs  (SWMEM  inputs)  within  it.  Since  his  routine  will  sup¬ 
plant  that  written  by  PREC0MP,  the  effect  of  any  input  expressions  coded  only 
on  *SWTCH  (*SWMEM)  cards  will  be  lost.  The  procedure  for  writing  SWINPT 
and  SWMEMN  is  almost  identical  and  is  outlined  in  the  following  steps: 

•  Place  the  card,  ^SWITCHES  n,  (*SWMEMCNT  n)  where  n  is  the 
number  of  SWTCH's  (SWMEM's)  being  used,  among  the  run-time 
data  cards. 

•  Write  the  first  four  cards  of  SUBR0UTINE  SWINPT 
(SUBR0UTINE  SWMEMN),  normally  written  by  the  pre¬ 
compiler,  exactly  as  they  appear  in  Appendix  D-3-b. 

•  Include  any  common  blocks,  dimension  statements,  function 
definitions,  and  computations  needed  to  define  the  switch  input 
expressions. 

•  Define  each  of  the  SWTCH  (SWMEM)  inputs  in  the  form 

VALUES(i)  =  input  expression^ 

•  Conclude  the  subroutines  by  writing  the  cards  RETURN  and  END. 

2.  User-Computed  SWTCH  and  SWMEM  Output 

The  switch  output  variables  SWTCH(i)  and  SWMEM(i„j)  are  passed  in 
the  common  block  SWTCHS  to  those  segments  of  the  coding  translated  by 
the  precompiler.  Therefore,  they  are  available  at  any  time  within  the 
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subroutines  DERIVS,  0UTPUT,  SWINPT,  and  SWMEMN.  SWCHi  and  SWMi, 
however,  are  computed  only  within  the  derivative  segment  and  appear  there 
only  if  *SWTCH  and  *SWMEM  statements  have  been  included,  [if  *SWTCH 
(*SWMEM)  cards  are  used  in  addition  to  user-written  SWINPT  (SWMEMN), 
PREC0MP  will  correctly  write  the  coding  into  DERIVS  to  define  SWCHi 
(SWMi),  even  though  the  input  expressions  will  be  lost.]  There  are  situations, 
therefore,  in  which  the  user  may  have  to  use  SWTCH(i)  and  SWMEM(i,  j)  to 
compute  SWCHi  and  SWMi  himself:  the  first  is  during  routine  usage  of 
*SWTCH  or  *SWMEM  when  the  user  wants  to  compute  SWCHi  or  SWMi  to  use 
outside  of  DERIVS,  and  the  second  is  when  he  has  written  his  own  switch 
input  routines,  uses  no  *SWTCH  or  *SWMEM  cards,  and  wants  to  use  these 
variables  anywhere. 

SWCHi  can  be  computed  in  any  routine  which  contains  the  common 
block  SWT  CHS  by  including  the  statements 

IF  (SWTCH(i)  .  GT.  0)  SWCHi  =  some  expression 

IF  (SWTCH(i)  .  LE.  0)  SWCHi  =  some  expression 

SWMi  can  also  be  computed  in  any  routine  containing  the  common  block 
SWTCHS  by  using  the  expression 


SWMi  =  SWMEM(i,  3)-SWMEM(i,  2)*(SWMEM(i,  1) -input) 


E.  HOW  THE  SWITCHES  WORK 

If  switches  of  any  form,  SWTCH,  SWMEM,  or  EVENTS,  are  used  in 
an  ESP  program,  their  inputs  are  processed  regularly  after  the  completion 
of  successful  integration  steps  to  see  if  any  switches  have  occurred.  To 
minimize  the  calculations  required  to  do  this,  the  inputs  are  all  defined  in 
the  separate  routines --SWINPT,  SWMEMN,  and  EVENTS--so  that  only 
these  routines  (and  not  DERIVS)  need  to  be  called  to  check  the  inputs.  To 
detect  switches  and  find  the  zero  crossings,  one  past  value  of  the  input  is 
always  saved. 
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Once  a  switching  has  been  detected,  a  modified  version  of  Wilkinson's 
method  (Ref.  2)  is  used  to  find  the  time  of  switching.  If  neither  the  saved 
value  nor  the  present  value  of  the  input  is  zero,  linear  interpolation  is  used 
four  successive  times,  testing  for  convergence  each  time.  If  convergence 
has  not  been  accomplished  after  four  iterations,  a  bisection  is  performed, 
and  the  linear  interpolation  is  repeated  in  sequences  of  four  plus  a  bisection 
until  convergence  is  achieved.  The  condition  for  convergence  is  that  the 
zero  crossing  be  found  within  an  interval  of  time  (see  the  explanation  of  the 
*HSW,  *HSWM,  and  *HSWE  cards.  Section  V-F). 

After  determination  of  the  first  switching  (if  more  than  one  occurred 
in  the  interval  T  -  H  to  T),  all  switches  which  would  switch  within  the 
accuracy  requirements  of  the  first  are  allowed  to  do  so;  then  the  solution  is 
restarted.  The  sequence  is  as  follows: 

1.  The  zero  crossing  is  found. 

2.  N0TIFY  is  called  if  EVENTS  is  used  and  any  have  occurred. 

3.  Printing  is  done  at  any  print  intervals  prior  to  the  switch  time. 

4.  Variables  to  be  used  in  restarting  the  integration  are  recomputed 
at  the  switch  time. 

5.  If  NDISPR  =  2,  the  output  data  prior  to  the  switching  (that  is, 
the  data  at  the  switch  time  but  before  the  effect  of  the  switch  has 
been  computed)  is  plotted  and  printed. 

6.  SWTCH(i)'s  and  SWMEM(i,4)‘s  are  set  to  the  proper  signs  and 
values  and  0.  5  is  added  to  those  that  have  switched. 

7.  DERIVS  is  called  to  evaluate  the  derivatives  with  the  0.5  flags 
on  the  switch  outputs  (see  Sections  V-A  and  V-B-2). 

8.  All  SWTCH  and  SWMEM  inputs  are  reevaluated  and  their  outputs 
updated  in  case  one  switch  has  toggled  another. 

9.  EVENTS  inputs  are  reevaluated  if  used. 

10.  The  0.5's  are  stripped  from  the  switch  outputs. 

11.  DERIVS  is  called  to  evaluate  the  derivatives  without  the  0.5's 
on  the  switch  outputs. 
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12.  Plot  data  is  stored  and  if  NDISPR  *  0,  print  occurs. 

13.  Integration  is  restarted  from  the  switch  time.  (See  flow  charts 
of  ESPCTL  in  Appendix  D-4-c-iii.  ). 

If  another  switch  occurred  later  in  the  same  integration  interval,  it 
will  be  detected  after  the  next  successful  integration  step,  and  the  above 
procedure  will  be  repeated. 

This  sequence  of  events  has  several  important  implications  for  pro¬ 
grams  in  which  two  or  more  SWTCH's  or  SWMEM's  occur  in  series: 

•  Two  SWTCH's  in  a  serits  (i.  e.,  the  first  triggers  the  second)  will 
produce  the  correct  output  for  the  second  SWTCH  in  the  second 
call  to  DERIVS.  Also,  SWTCH(i)  will  contain  an  accurate 

count  of  the  actual  number  of  switchings  to  date,  but  no  flag 
will  appear  on  the  second  SWTCH(i).  Thus  the  "IF  (SWTCH(i) 

.  EQ.  AINT (SWTCH(i))).  .  .  "  test  will  not  work  for  the  second 
switch. 

•  Two  SWMEM's  in  a  series  (i.  e.,  the  first  triggers  the  second)  will 
also  produce  the  correct  output  for  the  second  SWMEM  in  the 
second  call  to  DERIVS.  On  both  the  first  and  second  calls  to 
DERIVS,  the  integer  value  of  SWMEM (i,  4)  will  accurately  indi¬ 
cate  the  position  on  the  hysteresis.  However,  no  flag  will  ap¬ 
pear  on  the  second  SWMEM  and,  like  the  second  SWTCH(i),  the 
user  will  not  be  able  to  test  for  it. 

•  More  than  two  SWTCH's  or  SWMEM's  in  a  series  will  result  in 
the  first  two  being  detected  as  above,  and  the  next  two  being 
detected  and  reported  some  HSW  or  HSWM  interval  later,  and 
so  on. 

•  The  amount  and  accuracy  of  print  and  plot  data  in  the  neighbor¬ 
hood  of  switches  may  be  controlled  by  use  of  the  flag  NDISPR: 

NDISPR  *  0  Plot  data  is  stored  at  end  of  switch 

sequence,  but  no  print  occurs  at  the 
switch  time  unless  it  happens  to  coincide 
with  a  print  time. 

1  Plot  data  is  stored  and  print  occurs  at 
the  end  of  the  switch  sequence. 

2  Plot  data  and  print  data  are  stored  at  the 
switch  time  both  before  the  effect  of  the 
switch  is  calculated  and  afterward. 
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The  default  value  of  NDISPR  is  1.  To  change  this  the  user  must  include 


C0MM<jN/ NDISPR  /  NDISPR 


in  his  program  (preferably  in  ICC0MP)  and  set  NDISPR  to  the  appropriate 
value. 

F.  CONTROLLING  TIMING  ACCURACY  OF  DISCONTINUITIES 
(*HSW,  *HSWM,  AND  *HSWE) 

There  are  three  special  control  cards  for  controlling  the  allowable 
timing  error  in  determining  SWTCH's,  SWMEM's,  and  EVENT'S.  All  are 
optional  and  if  used  are  placed  among  the  run-time  data  cards,  that  is  after 
*ICCOMP. .  .  *ENDIC  and  before  *RUN,  in  any  order.  Their  formats  are 


*HSW 

hl 

h2 

...  h 

n 

$ 

*HSWM 

hl 

h2 

•  •  •  h 

n 

$ 

*HSWE 

hl 

h2 

...  hn 

$ 

where 


Default  values  of  all  h^  =  IE  -  6 

*HSW  (*HSWM,  *HSWE)  starts  in  column 


h.  can  be 
1 


a  constant  alone 
Hi  =  constant 
ALL  =  constant 


1 


(If  a  constant  appears  alone,  it  is  assumed  to  be  the 
control  for  the  next  SWTCH,  SWMEM,  or  EVENT.  ) 

$  is  a  required  terminator 
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EXAMPLES: 


*HSW 

produces : 

*HSWM 

produces : 


*HSWE 

produces : 


ALL=  1.  E  - 1 0  $ 

all  HSW(i)  =  1.  E-10 
HU1.E-2  H3  =  l.  E-10  $ 

HSWM(1  )=  1.  E  -2 
HSWM(2)=1.  E-6  [default] 
HSWM(3)=1.  E-10 
any  additional  HSWM(i)=l.E 

H2  =  l.  E-20  l.E-3  H6=l 


HSWE(2)=1.  E-20 

HSWE(3)=1.  E-3 

HSWE(6)=1.  E-10 

HSWE(  1  )=HSWE(4)  =  HSWE(5 ) 
HSWE(i)=l.  E-6  [default] 


-6  [default] 
.E-10  $ 


any  additional 
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SECTION  VI 


OUTPUT 


Output  from  ESP  may  take  a  variety  of  forms,  depending  upon  the  needs 
and  wishes  of  the  user.  Printed  output  may  be  automatically  formatted  by 
means  of  the  ESP  command  ’.-PRINT  or  it  may  be  tailored  to  the  user's  speci¬ 
fications  using  standard  FORTRAN.  Graphic  output  may  be  produced  on 
computer  printout,  on  microfilm  or  on  paper  by  the  Calcomp  pen  plotter 
(see  "IPD  Computing  Guide,  "  Ref.  4).  Data  files  may  be  written  onto  mag¬ 
netic  tapes  for  later  use  by  another  program  or  for  later  plotting.  Any  or 
all  of  these  modes  of  output  may  be  combined  in  any  one  program. 

Because  each  mode  of  output  has  advantages  for  particular  situations, 
the  following  sections  will  attempt  to  give  the  user  sufficient  information, 
not  only  to  easily  use  each  one,  but  also  to  help  him  decide  which  best  suits 
his  needs. 

A.  PRINTED  OUTPUT:  AUTOMATIC  FORMATTING 

The  fastest  and  easiest  way  for  the  user  to  obtain  printed  output  from 
his  program  is  with  the  SPRINT  command.  With  this  option,  the  user  needs 
only  to  name  the  variables  to  be  printed  and  the  labels  to  be  assigned;  he 
gives  no  print  formats.  Labels  and  values  are  printed  in  the  order  named, 
six  columns  to  a  page,  in  E-format.  All  labels  are  printed  and  then  all 
values,  so  that  if  there  are  more  than  six  of  each,  corresponding  labels  and 
values  will  not  be  adjacent,  although  their  correspondence  will  still  be  clear 
(see  example  below).  *PRINT  may  be  used  in  two  different  ways,  depending 
on  the  nature  of  the  values  to  be  printed. 

1 .  Printing  ESP  Variables  (*PRINT) 

If  all  the  values  to  be  printed  can  be  expressed  in  terms  of  T,  Y,  DY, 
PAR,  constants,  user  supplied  functions,  or  system  functions  such  as  SIN, 
COS,  and  SQRT  (i.e.,  no  other  variable  names  are  needed  to  define  the 
values  for  output),  then  the  fastest  and  easiest  way  to  obtain  printout  is  to  use 
the  ESP  command  card 
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SPRINT  label|  =  expression^  $. .  .  .  labeln  =  expression^  $  $ 


where 


SPRINT  starts  in  column  1 

label.  is  10  characters  (8  on  IBM)  or  less  (with  no  embedded 

blanks)  which  will  be  used  to  label  the  value  of  expres¬ 
sion  at  print  times 

expression.  is  any  F0RTRAN  expression  using  T,  Y,  DY,  PAR, 

constants,  and  system  functions  up  to  1206  charac¬ 
ters  (No  other  variables  may  be  used.  ) 

$  $  terminates  the  entire  *PRINT  statement  (No  continua¬ 

tion  marks  are  used  and  all  72  columns  may  be  used.  ) 


NOTE 

SPRINT  may  appear  only  once  in  a  program 
if  it  is  used  without  *0UTPUT  and  *END0UT. 


The  *PRINT  statement  will  be  translated  into  FORTRAN  statements 
which  will  be  written  as  part  of  SUBROUTINE  0UTPUT  in  the  form 

PRINT  (i)  =  expression^ 

The  labels  will  be  written  onto  a  ^HEADINGS  card  which  will  be  placed  on 
TAPE12  with  the  other  run-time  data  cards.  The  number  of  labels  found  on 
the  ^HEADINGS  card  at  execution  time  will  determine  the  number  of  values 
actually  printed  from  the  PRINT  vector  (see  SHEADINGS  card,  Section  VII-G-1.) 

EXAMPLE: 


SPRINT  TIME=T  $  SIGMA11=Y(1)*PAR(1)  $ 

SIGMA12  =  Y(2)*PAR(1)  $  SIGMA21  =  Y(3)*PAR(1)  $ 

SIGMA22=Y(4)*PAR(1)  $  PHI=Y  (5)*PAR(1)  $ 

PSI=Y(7)*PAR(1)  $  THETA  =  Y(6)*PAR(1)  $ 

OMEGA  1  =  Y(8)*PAR(1)  $  OMEGA2=Y(9)*PAR(l)  $ 

OMEGA  3=Y(10)*PAR(1)  $  THP=(Y(6)- Y(1  P)*PAR  (1)  $  $ 


produces  a  printout  that  looks  like  this: 


Tint 
PS  I 

1.  13j03..£*  J  3 

6.  45  7  0  ♦ u  3 

ji 

7,  L7Juli*z.*','j 

1*  30 J 3 uu£  +  j A 
7.672olu£Oj 


5I5MA11 
TnEI  A 

2.  05H34£*G1 
-5. 55 576 ?£ *01 

2.i05^78E*C'l 
-7 . 25  374 1 1»  G 1 

2.147414£*G1 

-7.856244E»G1 


S1GHA12 
OMEGA  1 

1 .24 15u2E*01 
-8  .4  16443E-03 

1.29QH8£*3t 
-8 .642948E-03 

1.3J5265£*0t 
-9.2163  35£-03 


SIGMA21 

0MEGA2 

9.999832E-39 
-5. 938054 £-32 

9.999832E-D9 

-5.935846E-02 

9.999832C-39 

-5.934609E-32 


SIGMA2  2 
0MEGA3 

9. 9998  22E-09 
4.5C9746f-03 

9.9998  J2F  -09 
4. 2545  0 1E-0  3 

9.999832E-C9 

3.992559E-C3 


PHI 
T  HP 

-1.59  75o8E*0  3 
-6.655762E*01 

-  1 . 74  57  39E ♦ 3  0 
-7.255741E*91 

1.89C354£»33 

-  7 , 856244E*  0 1 


2.  Printing  User  Variables  or  Computing  Output 
(^OUTPUT.  .  .  *END0UT) 

To  print  the  values  of  expressions  which  contain  variable  names  other 
than  T,  Y,  DY,  PAR,  system  functions  or  function  subprograms,  the  ESP 
command  cards,  *0UTPUT  and  *END0UT  must  be  used  to  signal  the  beginning 
and  the  ending  of  the  output  computations  and  printing  commands.  All  state¬ 
ments  between  these  two  cards  will  be  written  into  the  output  subroutine,  and 
the  user  will  find  it  helpful  to  remember  that  this  program  segment  is  a 
separate  subroutine  as  he  decides  what  he  may  and  may  not  do  within  it.  In 
general,  between  ^OUTPUT  and  *END0UT,  any  FORTRAN  or  WHELP  may  be 
used  to  define  user  output  variables,  and  *PRINT,  as  defined  above,  may  be 
used  to  print  them;  but  certain  rules  must  be  followed: 

•  If  *PRINT  and  *0UTPUT.  . .  *END0UT  are  both  used,  *PRINT 
may  only  be  used  between  *0UTPUT  and  *END)t)UT,  but  it  may 
appear  as  many  times  as  necessary  here  as  long  as  the  total 
number  of  print  variables  specified  does  not  exceed  60. 

•  User  variables  appearing  in  *PRINT  expressions  must  be  defined 
before  SPRINT. 

•  FORTRAN  rules  regarding  the  sequence  of  declarative  and 
executable  statements  within  the  section  must  be  followed. 
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The  general  format  of  this  section  is 


Ik 


OUTPUT 

FORTRAN  and  WHELP  statements  (in  any  order  consistent  with 

FORTRAN  rules) 

*PRINT  label,  =  expression,  $...  label  =  expression  $  $ 

11  n  n 

FORTRAN  and  WHELP  statements 
*END0UT 

EXAMPLE: 

*0UTPUT 

C0MM0N/BL0CKA/A,  B,  R,  0MEGA 
'-TDECLARE  A (6,  3)  B{3 . 6)  R(6)  X(6)  $  [WHELP  statement] 

CALL  DIDDLE  (T,  Y,  THETA ,  A  LPHA ) 

X  =  A  *  B  *  R/  (THETA)  $  [WHELP  statement] 

Z  =  THETA  +  T  *  A  LPHA 
♦PRINT  TIME  =  T  $  X1=X(1)  $  X2=X(2)  $ 

X3  =X  (3 )  $  THETA  =  THE  TA  $  Z  =  Z  $  $ 

(additional  F0RT RAN) 

*END0UT 

Notice  that  the  subroutine  DIDDLE  must  be  supplied  by  the  user, 
and  that  the  variables  A,  B,  R,  and  0MEGA  must  be  defined  elsewhere. 

3 .  Accuracy  of  Printed  Values 

Since  the  print  interval  and  stepsize  are  unrelated,  it  is  generally 
necessary  to  evaluate  the  solution  just  for  printout.  Adams  integration  uses 
the  code  and  method  outlined  in  Ref.  6  ,  and  the  solution  (Y's)  and  deriva¬ 
tives  (DY's)  are  interpolated  from  the  difference  table  kept  internally.  For 
predictor-corrector  and  Runge-Kutta  integration,  in  order  to  avoid  calling  the 
integration  routine  just  for  printout,  Hermite  interpolation  (Ref.  3),  which  uses 
the  functional  values  and  their  derivatives  at  three  points,  is  employed  to 


evaluate  the  solution  (Y's)  and  derivatives  (DY's)  at  intermediate  points. 

Thus,  the  proper  T  (time)  and  its  corresponding  Y's  and  DY's  are  automati¬ 
cally  passed  to  0UTPUT  at  a  print  time,  and  these  values  and  any  other  output 
variables  computed  using  only  these  and  constants  will  be  correct  for  the  time 
g  i  ve  n . 

However,  since  in  the  process  of  integration  ESP  generally  oversteps 
the  print  time  and  "backs  up"  to  print,  PAR's  and  time  dependent  variables 
which  may  have  been  passed  to  0UTPUT  by  a  user- supplied  common  block 
may  not  correspond  to  the  print  time.  The  best  way  to  avoid  this  problem  is 
to  recompute  these  variables  within  :':0UTPUT.  .  .  *END0UT,  so  that  they  will 
always  reflect  the  actual  values  at  the  print  time.  If  this  involves  much 
computation,  the  user  may  want  to  test  PRINT ( 1 ) ,  as  shown  in  the  example 
in  Section  VI-C,  to  ensure  that  these  values  are  recomputed  only  for  printout. 
(Since  plot  data  storage  occurs  at  the  end  of  each  successful  integration  step, 
or  pair  of  steps,  irrespective  of  the  print  interval,  data  stored  for  plotting 
will  be  consistent  and  the  user  need  not  concern  himself  with  this  problem.  ) 

B.  GRAPHIC  OUTPUT 

To  produce  graphic  (plotted)  output  from  an  ESP  run,  the  user  has  two 
tasks:  storing  the  data  to  be  plotted  and  specifying  how  it  is  to  be  plotted. 

The  simplest  way  to  do  these  tasks  is  to  store  data  into  the  vector  PE0T  and 
then  use  the  *GRAPH  command  (explained  below)  to  plot  it.  An  alternative 
way  is  to  write  all  plot  data  onto  a  magnetic  tape  or  disk  file  and  then  plot  it 
using  some  other  plotting  routine  (see  "IPD  Computing  Guide,  "  Ref.  4).  The 
first  method  is  the  simplest  and  most  satisfactory  for  most  user  needs.  The 
latter  is  useful  mainly  for  very  time-consuming  program,  runs  where  it  is 
desirable  to  have  output  data  available  for  repeated  plotting  or  study  without 
the  necessity  of  rerunning  the  program.  (Refer  to  Section  VI-D  on  Data  File 
Output.  )  It  is  also  possible  to  produce  overlays  using  data  generated  in  dif¬ 
ferent  runs  (see  Appendix  F  on  Multiple  Runs.  ). 


-■  •--■ftwtarii  i  • 
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1. 


Storing  Plot  Data 

All  data  to  be  plotted  by  *GRAPH  must  be  stored  by  the  user  in  the  vec¬ 
tor  PL0T.  The  simplest  way  to  do  this  is  within  the  *PRINT  statement  used 
with  or  without  *0UTPUT.  .  .  *END0UT,  as  follows: 

*  PRINT  TIME=PL0T(1)  =  T  $  RATE=PL0T  (2)  =  Y(1)  $ 

0MEGA  =  PL0T  (3  )  =  Y  ( 1 )  +  Y(2)  $  $ 

This  statement  accomplishes  both  printing  and  labeling  as  explained 
above  and  storing  of  data  for  plotting.  It  generates  the  following  statements 
in  the  FORTRAN  version  of  SUBROUTINE  OUTPUT. 

PRINT  ( 1  )=PL0T  ( 1  )  =  T 
PRINT  (2  )=PL0T (2)  =  Y(  1 ) 

PRINT  (3  )  =  P  L0T  {3  )  =  Y  (1 )  +  Y(2) 


The  user  may  also  store  the  plot  data  himself  if  he  is  not  using  *PRINT 
or  wishes  not  to  bother  adding  to  or  changing  his  *PRINT  cards.  To  do  so 
he  simply  adds  the  FORTRAN  statements  needed  to  the  *0UTPUT .  .  *END0UT 
section  of  his  program.  The  example  above  would  then  become 
*0UTPUT 

*PRINT  TIME=T  $  RATE  =  Y(1)  $  0MEGA=Y(1)  +  Y(2)  $  $ 

PD0T(1)=T 
PL0T(2)  =  Y(1) 

P  LOT  (3  )  =  Y  (1 )  +  Y (2 ) 

*END0UT 


WARNING 

Do  not  attempt  to  store  plot  data  in  this  way  in  other  parts  of 
the  program  (anywhere  outside  of  *0UTPUT.  .  *END0UT) ;  array 
PLOT  will  not  be  recognized  and  the  timing  of  data  storage  will 
be  wrong. 


The  assumed  or  default  number  of  plot  variables  which  may  be  stored  in 
this  way  is  10.  If  more  than  10  variables  are  to  be  stored,  a  *MAXPL0TS 
card  must  be  placed  among  the  run-time  data  cards,  that  is,  after  *ICC0MP.  .  . 
*ENDIC,  but  before  *RUN.  The  format  is 
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*MAXPL0TS  n 


where 

*MAXPL0TS  starts  in  column  1 

n  is  an  integer  (<  100)  specifying  the  maximum  number  of  plot 

variables  to  be  used 


During  execution  of  the  program,  ESP  uses  a  file  (TAPE11)  to  store 
the  data  placed  in  the  PL0T  vector,  storing  a  point  at  the  end  of  each  suc¬ 
cessful  integration  step  or  pair  of  steps  and  at  each  switch  time.  At  plotting 
time  ESP  selects  a  representative  sample  of  this  data  for  actual  plotting  (up 
to  1000  values  per  variable),  and  in  most  cases  no  significant  loss  of  infor¬ 
mation  occurs.  However,  if  the  user  wishes  to  have  greater  control  over 
the  number  of  points  or  the  intervals  over  which  they  are  plotted,  PR0GRAM 
ESPPL$T  may  be  used  instead.  (See  Section  VI-D-1.) 

2.  Plotting  Output  (*GRAPH) 

Printer,  film,  and/or  pen  plots  are  obtained  by  placing  *GRAPH  com¬ 
mands  after  the  *RUN  which  computes  the  results  to  be  plotted.^  The 
general  form  is 


where 


nx  (1  $  n^  £  100)  is  the  PL0T  subscript  of  the  desired  x  variable 


n^.  (1  ^  n^  ^  100)  is  the  PL0T  subscript  of  the  desired  y  variable 
([size],  [grid],  [scaling],  and  [type]  are  optional,  and  may  appear  in 


any  order  after  n  and  n 
*GRAPH. )  X 


but  must  appear  on  the  same  card  with 


Plotting  is  accomplished  by  the  subroutine  GRAPH;  for  more  complex  needs 
and  some  further  options,  the  user  is  referred  to  the  writeup  of  this  sub¬ 
routine. 
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More  specifically 


where 


[size] 


’SMALL 
SIZExxyy 
OVERLAY 
.  OVER  LAY  1. 


Default  choice  is  SMALL 


SMALL 


SIZExxyy 


OVERLAY 


OVERLAYl 


implies  a  6  X  10-in.  printer  plot  or  a  10  X  15  in. 
hardcopy  or  film  plot 

implies  an  xxXyy-in.  plot  (10  X  10-in.  or  10  X  15-in. 
nicely  fills  a  linear  microfilm  grid;  this  will  be  turned 
sideways  and  possibly  cover  more  than  one  page  on 
the  printer  plots. ) 

implies  size  and  scaling  as  on  the  previous  graph.  A 
new  plot  will  result  on  the  printer  while  a  true  over¬ 
lay  will  be  made  on  pen  or  film  plots  (See  TYPES). 

implies  the  same  as  OVERLAY  except  that  the  data 
is  completely  rescaled  and  on  pen  and  film  a  new  y- 
axis  scale  is  placed  at  the  right  end  of  the  graph. 


[grid]  =  [GRIDggg] 


whe  re 


Default  choice  is  GRID3A1 


ggg  is  the  three-character  number  specifying  the  type  of  plot  grid 
to  be  used  (see  "IPD  Computing  Guide,  11  Ref.  4) 

Printer  plots  are  always  linear-linear 

If  grid  is  semilog  or  log-log,  film/pen  plots  are  made  using  log,n  of 
the  X  or  Y  data  or  both 


6-8 


[scaling] 


SCALE 


AUTO 

X  XDEL 
o 


AUTO 

Y  YDEL 
o 


where 


Default  choice  is  AUTO  for  both  axes  (automatic  scaling  based  on 
actual  data  stored) 

If  SCALE  appears,  the  parameters  within  must  be  specified  for 
both  X  and  Y  11 


AUTO  produces  automatic  optimized  scaling  based  on  actual  data  for 
X-axis  (Y-axis) 

X  (Y  )  is  the  minimum  scale  value  to  be  used  for  X  data  (Y  data) 
o  o 

XDEL  (YDEL)  is  the  absolute  value  of  the  difference  between  one 
scale  annotation  and  the  next.  For  a  10  X  10-in.  plot,  XDEL  = 

(XMAX  -  X0)/10. 

[type]  =  [TYPEt1[t2][t3] 

whe  re 

Default  choice  is  TYPEP 

P  for  printer  plots 
or 

S  for  printer  plots  to  be  overlaid 
F  for  film  plots 
C  for  Calcomp  pen  plots^ 

Any  combination  of  P,  F,  and  C  may  be  used  (Example: 

TYPEFC,  TYPEP,  TYPECFP) 

S  may  be  used  in  place  of  P  if  Printer  plot  overlays  are  to  be  made. 

If  TYPES  is  specified  on  the  *GRAPH  card,  no  printer  plot  will  be 
produced  for  that  card,  but  the  plot  data  will  be  stored.  Thus, 
utilizing  TYPES  on  all  but  the  final  overlay  card  will  produce  a 
single  printer  plot  containing  all  overlays  to  be  produced. 

See  "IPD  Computing  Guide,  "  Ref.  4,  for  other  steps  required  to  obtain  pen 
plots . 


t.  = 

i 
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EXAMPLE: 

The  following  sequence  of  *GRAPH  cards  will  produce  a  single 
printer  plot  with  three  graphs  on  it. 


-GRAPH  1 
-GRAPH  1 


TYPES 

TYPES  OVERLAY 


*GRAPH  1  4  OVERLAY 


[title (s )] 


plot  title 
X  -title 
Y-title 


whe  re 


Default  choice  is  blanks 

Each  title  is  any  character  string  appearing  in  columns  1  through  50. 
(Column  1  should  be  blank  to  avoid  possible  misidentification.  ) 

All  titles  are  optional  but  must  appear  in  the  order  given,  each  on  a 
separate  card  (To  delete  one,  substitute  a  blank  card.  ) 


EXAMPLES: 


1.  -RUN 


-GRAPH  1  4 

This  will  produce  a  printer  plot  of  the  Y-axis  data  stored  in  PL0T(4) 
versus  the  X-axis  data  stored  in  PL@T(1).  It  will  make  a  SMALL  (6  X  10-in.  ) 
plot  with  automatic  scaling  and  no  titles.  This  represents  the  minimum 
specifications  to  produce  a  plot. 


2.  -RUN  7 
-GRAPH 


l  5  SIZE1510  TYPEPF  SCALE  AUTO  5.0  50 

ERROR  IN  ANGLE  OF  ATTACK 
DESIRED  ANGLE 

DIFFERENCE  BETWEEN  GOA  L  AND  ACTUAL 
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This  example  shows  considerable  user  control  of  plot  parameters. 

It  will  produce  a  linear  plot  on  printer  and  film  with  an  X-axis  15-in. 

long  and  automatically  scaled.  The  Y-axis  will  be  10  in.  with  Yq  at 

5.0  and  Y  at  505.  0.  "ERROR  IN  ANGLE  OF  ATTACK11  will  be 
max 

printed  across  the  top  of  the  plot,  "DESIRED  ANGLE"  along  the  X-axis 
and  "DIFFERENCE..."  along  the  Y-axis. 

C.  PRINTED  OUTPUT:  USER  FORMATTED 

At  any  time  the  user  wishes  to  produce  output  printed  to  his  own  speci¬ 
fications,  he  may  do  so  by  simply  adding  FORTRAN  print  (or  write^) 
statements  and  their  corresponding  formats  to  any  section  of  his  program, 
even  if  he  is  also  using  *PRINT  or  *0UTPUT.  .  .  *END0UT.  Since  con¬ 
siderable  care  must  be  exercised  in  controlling  the  frequency  of  printout 
produced  in  this  manner,  it  is  suggested  that  *PRINT  be  used  as  much  as 
possible  during  program  development  and  debugging  stages  because  of  its 
speed  and  ease,  and  that  F0RTRAN  print  statements  be  added  only  when  a 
more  specific  format  is  needed  for  production  runs.  In  adding  print  state¬ 
ments  the  user  should  consider  the  following  points: 

1.  Since  each  segment  of  the  user's  coding  results  in  a  separate 
subroutine,  the  user  must  be  careful  to  print  in  each  segment 
only  those  values  known  within  the  section;  that  is,  those 
values  passed  to  it  by  ESP,  defined  within  the  segment,  or 
passed  to  the  segment  by  user  -  supplied  common  block. 

2.  A  print  statement  placed  within  *ICC0MP.  .  *ENDIC  will  be 
executed  once  at  the  beginning  of  each  run,  since  that  is  the 
only  time  SUBROUTINE  ICC0MP  is  called.  This  is  the  logical 
place,  therefore,  to  print  variables  which  are  constant  for  the 
particular  run.  (Refer  to  Section  II-B-STEP5  for  a  list  of  the 
inputs  and  constants  ESP  prints  automatically  at  program  ini¬ 
tiation.  ) 


€ 

Since  the  main  program  written  by  PREC0MP  does  not  assign  a  TAPEn  type 
name  to  the  Output  file,  if  WRITE  statements  are  used,  the  user  must  write 
his  own  main  program  in  order  to  set  TAPEn  =  OUTPUT. 


Placed  within  *DERIVS.  .  .  *ENDDERIVS,  print  statements  will  be 
executed  every  time  the  derivatives  are  evaluated,  which  is 
several  times  per  step  at  possibly  decreasing  times,  the  exact 
number  depending  on  the  integration  algorithm  used.  Printout 
of  this  type  may  be  useful  for  specific  debugging  information, 
but  will  be  unwieldy  and  confusing  for  general  output. 

Data  which  is  to  be  printed  only  at  the  completion  of  a  run  may 
be  printed  from  the  MAIN  program.  The  user  would  need  to 
(1)  place  a  ^RETURN  card  at  the  end  of  the  run-time  data  cards, 
to  cause  program  control  to  return  to  MAIN;  (2)  write  his  own 
MAIN  (see  APPENDIX  D-3-b)  adding  common  blocks  to  transmit 
the  data  to  be  printed;  and  (3)  add  print  coding  to  MAIN  between 
CALL  ESPII  and  END. 

Within  *0UTPUT.  .  *END0UT  is  generally  the  best  place  to  add 
print  statements,  since  the  output  subroutine  is  called  at  regular 
and  predictable  intervals,  namely,  once  at  the  end  of  each  inte¬ 
gration  step  for  plotting  purposes  and  once  at  each  print  time  for 
printing  purposes.  As  long  as  printing  is  accomplished  by 
*PRINT  commands,  ESP  checks  to  see  that  printing  occurs  only 
at  the  specified  print  intervals.  When  the  user  writes  his  own 
print  statements  and  formats,  however,  he  must  be  the  one  to 
see  that  they  are  executed  at  the  proper  times.  To  do  this  he 
simply  tests  PRINT(l)  (before  setting  it!).  PRINT(l)  will  be 
4HPL0T  if  OUTPUT  is  being  called  for  plotting  purposes  or 
5HPRINT  if  it  is  called  at  a  print  time. 


EXAMPLE: 

^OUTPUT 

If  (PRINT (1 )  .  NE.  5HPRINT)  GO  TO  5 
PRINT  100,  T,  Y(l),  DY (1 ) 

100  FORMAT  (1H0,  3E20.  8) 

5  CONTINUE 


*  PRINT  . 


*END0UT 
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D. 


DATA  FILE  OUTPUT 


In  addition  to  printed  and  plotted  output  from  an  ESP  program,  the 
user  may  wish  to  output  his  data  on  a  magnetic  tape  that  he  can  reuse  after 
the  termination  of  his  job.  Two  common  reasons  for  needing  this  capability 
are  the  desire  to  use  the  output  as  input  to  another  computer  program  and  the 
wish  to  ensure  that  data  from  a  particularly  time-consuming  run  is  saved  for 
plotting.  In  either  case,  there  are  two  alternative  ways  to  save  these  results, 
either  by  saving  TAPE11,  the  file  on  which  ESP  automatically  writes  all  plot 
output,  or  by  saving  a  file  created  and  named  entirely  by  the  user.  TAPE11 
is  easier  to  produce,  but  because  the  data  on  it  is  both  packed  and  blocked, 
a  special  program,  ESPPL0T,  must  be  used  to  plot  it  (see  below).  On  the 
other  hand,  a  user-created  file  can  be  tailored  specifically  for  its  later  use, 
but  requires  more  effort  to  generate. 

1 .  Data  Written  onto  TAPE  11 

ESP  automatically  writes  ail  data  stored  in  PL0T  (either  via  SPRINT 
statements  or  by  PL($T(i)  =  expression,  refer  to  Section  VI-B-I)  onto  a 
logical  file  named  TAPE11,  which  is  saved  until  run  completion,  plotted  from 
by  *GRAPH,  and  then  released.  To  save  this  data  after  job  termination, 
the  user  must  transfer  it  to  a  magnetic  tape.  Once  stored  on  magnetic  tape, 
the  data  may  be  plotted  at  any  time  using  program  ESPPLX^T. 

ESPPL0T  is  a  compiled  (binary)  main  program  which  is  to  be  run  at 
some  time  later  than  the  ESP  run  whose  data  it  is  to  plot.  It  may  be  used  to 
plot  up  to  10,000  points  and  allows  the  user  to  determine  the  time  intervals 
to  be  plotted.  Every  point  or  every  nth  point  may  be  plotted,  and  symbols 
may  be  placed  on  the  plots  at  any  At  interval  specified  by  the  user.  Steps 
required  for  its  use  are  the  following: 

a.  Make  certain  that  the  first  element  of  the  array  PL0T 
stored  by  the  ESP  run  is  time  (T). 

b.  Request  the  actual  tape  containing  the  data  and  assign  it 
the  logical  name  TAPE  11. 


c.  Execute  program  ESPPLOT  (See  control  cards  in  example  below.) 

d.  Write  the  necessary  -'-Control  cards  according  to  the  formats 
given  below. 

ESPPLOT  -‘CONTROL  CARD  FORMATS:  All  '-'‘Control  cards  must  start  in 
column  1.  Non-*  cards  may  occupy  columns  2  through  72,  except  title  cards, 
which  are  limited  to  columns  1  through  50. 


*MAXPLOTS  n 


*N LOCAL n 


TMAX  m  T 


(n  must  be  the  same  as  n  on  the 
*MAXPLOTS  card  appearing  in  the 
ESP  run.  Default  value  is  10.  ) 


(Every  nth  data  point  is  to  be  plotted. 
Default  value  is  1.  May  be  changed  at 
any  time .  ) 


-  .  .  T 

max,  '  * '  min 
1  m 


=  IMAX  m  ALL  T  .  T 

min  max 


*TICKPLOTS  At 


(The  next  m  sets  of  -‘GRAPHS  are  to  be 
plotted  for  the  m  sets  of  intervals  given. 
^minn  and  ^max  apply  to  the  nth -‘GRAPH 
encountered.  If  ftie  m  sets  of  graphs  are 
all  to  be  plotted  for  the  same  Tmin  and 
Tmax’  enter  ALL  followed  by  Tmin  and 
Tmax.  To  change  this,  use  ^RETURN 
first:  see  example.) 


(The  symbol  V  will  be  placed  on  the 
plotted  line  every  At  seconds.  The 
default  is  no  symbol.  ) 


*GRAPH  nx  ny  [type]  [size]  [grid]  [scaling] 
Title] 
x-title] 
y -title ] 


(See  *GRAPH  descrip¬ 
tion,  Section  VI-B-2 
above.  ) 
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:  RETURN 


(End  this  segment  of  instructions  and 
begin  another. 

*IMAX,  *TICKPIjOTS,  and  *GRAPH 
must  be  reentered. 

*N LOCAL  will  retain  its  value  but  ma 
be  changed. 

*MAXPLOTS  will  retain  its  value.  ) 


EXAMPLE: 


ATTACH,  ESPPLOT ,  2ESPPLOT. 
ATTACH,  SUBLIB,  2ESPFTN. 

ATTACH,  PLOTLB,  3FTNPLOT LIB. 
LIBRARY,  SUBLIB,  PLOTLB. 

(Request  TAPE  11  saved  by  ESP  run) 
ESPPLOT. 

HARDCPY. 

EOR 

*MAXPLOTS  12 
*NLOCAL  10 
*TICKP  LOTS  0.2 
*IMAX  2  ALL  0.  10. 

*GRA PH  1  9  TYPEF 

ETA  VS  TIME 
“•'GRAPH  9  10 

ETADOT  VS  ETA 
^RETURN 
*NLOCA  L  20 

-IMAX  3  50.  100.  50.  100.  0. 

*GRAPH  1  2  TYPEF 

*TICKPLOTS  5. 

*GRAPH  1  3  TYPEF  OVERLAY 

THETA  AND  TAU  VS  TIME 
*NLOCAL  10 
*GRAPH  2  5  TYPEF 

X  VS  THETA 
THETA 
X 

*RETURN 

EOF 


s  _ 


In  the  above  example  ESPPLOT  expects  TAPE  II  to  have  12  variables 
stored  on  it,  of  which  the  first  is  time.  Selecting  every  10th  data  point  and 
placing  a  tick  every  0.2  seconds,  it  will  produce  filmplots  of  variables  9  vs 
time  and  10  vs  9  from  T  =  0.  to  T  =  10.  seconds  and  label  them  ETA  VS  TIME 
and  ETADOT  VS  ETA,  respectively. 

A  ^RETURN  card  appears  next  so  that  IMAX  can  be  changed.  Note 
that  IMAX  in  this  section  specifies  3  time  intervals  and  therefore  3  -''-GRAPH 
cards  follow  it.  However,  only  2  plots  will  result,  as  the  first  two  are 
overlaid.  Selecting  every  20th  plot  point,  it  will  plot  variable  2  vs  time  from 
time  =  50.  to  100.  with  no  tick  marks.  Then  it  will  overlay  variable  3  vs 
time  on  the  same  grid  with  tick  marks  every  5.  seconds  and  label  this  plot 
THETA  AND  TAU  VS  TIME.  Finally,  using  every  10th  data  point,  it  will 
plot  variables  5  vs  2  from  time  =  0.  to  100.  with  tickmarks  every  5.  seconds, 
and  this  plot  will  have  X  VS  THETA  on  top,  THETA  on  the  X-axis  and  X  on 
the  y-axis.  A  final  ^RETURN  card  indicates  the  end  of  the  job. 

2.  Data  Written  onto  User-Named  File 

To  produce  a  saved  data  file  which  matches  some  particular  user  for¬ 
mat,  the  user  must  declare  and  write  his  own  file.  To  do  this  the  following 
steps  are  necessary: 

•  Request  that  a  magnetic  tape  be  allocated  to  the  ESP  program 
and  given  a  logical  file  name  by  using  the  proper  control  cards. 

•  Declare  the  logical  file  name  by  supplying  the  main  program 
(ESP's  PREC0MP  normally  writes  it),  and  adding  the  logical  file 
name  to  the  PROGRAM  card  (see  example  below), 

•  Write  desired  data  on  the  file  by  using  an  unformatted  WRITE(lfn) 
list  statement  within  ^OUTPUT.  .  .  *ENDOUT.  Remember  that  the 
output  subroutine  is  called  at  various  times  for  plotting  and 
printing  (refer  to  Section  VI-C),  and  that  any  time  the  user  does 
his  own  WRITE  statements  he  must  consider  the  frequency  with 
which  he  wants  to  "write"  and  allow  for  it  by  testing  PRINT  (1) 
and  acting  accordingly. 
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EXAMPLE: 

[control  cards] 

[control  card  additions  to  request  and  save  a  tape  and  to  name  it 
TAPE20,  according  to  current  operating  system  manual.] 

7-8-9 

PROGRAM  MAIN  (TAPE  12,  TAPE1 1, 0UTPUT,  INPUT  =TAPE  12,  TAPE20) 
EXTERNAL  DERIVS,  ADAMS,  ADMNTP 
CALLESPII  (DERIVS,  ADAMS,  ADMNTP) 

END 

[User  supplied  subroutines] 

*  DERIVS 

[Equations  defining  derivatives  and  switches] 

*ENDDERIVS 

Output 

[Equations  defining  outputs] 

If  (PRINT  (1)  .  NE.  5HPRINT)  GO  TO  5 
WRITE (20)  T,  Y(l),  Y(2),  Y(3),  0MEGA,  THETA 
5  CONTINUE 

[More  equations  defining  output,  *PRINT,  etc.] 

*END0UT 


Notice  that  the  MAIN  program,  here  supplied  by  the  user,  is  identical 
to  the  one  normally  written  by  PREC0MP  (see  Appendix  D-3-b)  except  for 
the  addition  of  TAPE20  to  the  PR0GRAM  card,  and  that  the  MAIN  program 
precedes  all  user  subroutines  and  the  *DERIVS  section.  Also  notice  that 
the  WRITE  statement  is  placed  inside  the  output  routine  (Do  not  place  it  in 
the  derivatives  section!)  and,  in  this  case,  that  it  will  be  executed  only 
when  0UTPUT  is  called  for  print  purposes,  which  will  be  at  the  print  inter¬ 
vals  specified  on  the  *RUN  card  and  at  switch  times  if  NDISPR  >  0.  This 
example  will  produce  a  file,  TAPE20,  with  six  values  per  record,  one 
record  for  each  print  time,  which  will  exist  on  magnetic  tape  after  run 
termination  and  which  may  be  read  into  another  program  for  input  or  plot¬ 
ting  by  an  unformatted  READ  statement. 
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SECTION  VII 


INPUTS 


The  inputs  needed  for  an  ESP  program  may  include  any  or  all  of  the 
following:  (1)  the  number  of  derivatives,  (2)  the  starting  and  ending  values 
for  the  independent  variable  T,  (3)  the  print  intervals )  desired,  (4)  the  initial 
conditions  (initial  values)  for  the  Y's,  (5)  variables  which  determine  solution 
accuracy,  (6)  miscellaneous  inputs  which  influence  program  control  and 
format,  and  (7)  any  user  variables  which  may  be  needed  to  compute  initial 
conditions,  derivatives,  or  switching  characteristics.  The  number  of  inputs 
required  and  the  manner  in  which  they  are  supplied  will  depend  upon  the  com¬ 
plexity  of  the  user's  problem  and  the  degree  of  flexibility  he  wishes  to  build 
into  the  program. 

A.  NUMBER  OF  DERIVATIVES.  START/STOP  TIMES,  AND 

PRINT  INTERVALS  (*RUN) 

The  only  inputs  required  by  every  ESP  job  are  the  exact  number  of 
derivatives,  the  initial  and  final  values  of  the  independent  variable  and  the 
printing  interval(s).  These  are  specified  on  the  *RUN  card,  which  is  re¬ 
quired  for  every  run  and  is  always  the  last  card  except  for  *GRAPH,  *ST($P 
or  ^RETURN  (refer  to  Section  VII-G-3).  Notice  that  the  print  interval  may 
be  changed  several  times  during  the  run,  and  that  t£  >  t  .  (Refer  to 
Appendix  F-l-e  for  directions  on  running  the  solution  backward.)  The 
format  is 

*RUN  neq  t  hprint.  t  hprint.,  ....  $ 

_ ° _ 1  *1 _ L  *2 _ 

where 

neq  is  an  integer  constant  specifying  the  number  of  dependent 

variables  (Y^).  It  should  equal  the  largest  subscript  of  DY 
and  must  be  corrected  if  derivatives  are  added  or  deleted. 

t  is  a  constant  specifying  the  initial  value  of  the  independent 

J  variable  (T)  (typically  0.) 


h 

P 


rint. 

1 


t 


f. 


1 


is  the  printing  interval  to  be  used  untiL  the  independent 
variable  (T)  reaches  t^ 

i 

is  either  the  value  of  T  at  which  the  solution  is  to  stop  or, 
if  it  is  followed  by  hprintn  tf  ,  it  is  the  time  at  which  the 
print  intervaL  is  to  be  changed. 


$  is  optional  and  terminates  the  information  on  the  card 


EXAMPLES: 

1.  -RUN  5  0.  0.5  10.  $ 

This  runs  the  solution  from  0.  to  10.0,  printing  every  0.5 
seconds,  solving  for  ’  dependent  variables  (Y. 's). 

2.  -RUN  12  2.0  0.25  8.0  0.5  20.0  $ 

This  solves  for  12  dependent  variables,  starting  at  2.0  seconds, 
and  printing  every  0.  25  second  until  T  reaches  8.  0  and  then 
printing  every  0.  5  second  until  T  reaches  20.0  seconds,  at 
which  time  the  solution  stops. 


B.  INITIAL  CONDITIONS:  KNOWN  CONSTANTS  (*1V) 

If  the  user  wishes  the  initial  values  of  all  the  dependent  variables  (YUs) 
to  be  zero,  he  does  nothing.  If  he  wishes  any  to  be  nonzero  constants,  the 
simplest  way  to  input  them  is  on  the  '-IV  card,  which  is  placed  among  the 
run  time  data  cards.  IV’s  retain  their  value  until  they  are  reset  on  another 
*IV  card  or  within  the  user's  coding.  The  format  is 

*IV  ivl  iv2  -  ivn  $ 


whe  re 

iv.  is  either  a  constant  alone  or  Yi  =  constant.  Values  are 

separated  by  blanks.  If  a  constant  appears  alone,  it  is 
assumed  to  be  the  initial  condition  of  the  next  dependent 
variable.  If  no  value  is  given  for  a  dependent  variable, 
it  is  assumed  to  be  zero. 

$  is  a  required  terminator 
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EXAMPLE: 

*IV  0.  Y6=3 . 0  6.  5  8.2  Y10=3.0  $ 

produces:  Y(1)=Y(2)=Y(3)  =  Y(4)=Y(5)=0.  0, 

Y(6 )=3 . 0,  Y (7)=6 . 5,  Y(8)=8.2, 

Y(9)=0.  0,  Y(10)=3.  0  and  any  additional 
Y(i)  =  0.  0. 


C.  USER  PARAMETERS:  KNOWN  CONSTANTS  (*PAR) 

An  array  named  PAR  of  length  100  is  created  by  ESP  and  automatically 
passed  by  the  statement  C0MM0N/PARS/PAR(1OO)  to  the  derivative,  output, 
switching,  and  initial  computations  routines  for  the  convenience  of  the  user. 
Its  main  purpose  is  to  permit  the  user  to  introduce  parameters  into  his  equa¬ 
tions  without  having  to  worry  about  how  they  will  be  defined  or  passed  to 
each  section,  and  to  permit  their  values  to  be  easily  changed  from  run  to  run. 

PAR  values  may  be  computed  in  ICC0MP  (see  below),  but  if  known  are 
most  easily  input  to  the  PAR  array  on  a  *PAR  card  which,  like  the  *IV  card, 
is  placed  at  the  end  of  the  program  but  before  the  *RUN  to  which  it  applies. 
PAR's  retain  their  values  until  they  are  reset  on  another  *PAR  card  or  within 
the  user's  coding.  They  may  be  used  at  any  time  by  name,  i.e.,  PAR(3). 


WARNING 

PAR's  must  not  be  functions  of  T.  Nonconstant  PAR's  used  as 
inputs  to  SWTCH's  or  SWMEM's  will  cause  errors  in  the  deter¬ 
mination  of  switching  times,  and  use  of  PAR  to  pass  nonconstant 
values  to  be  printed  will  cause  discrepancies  in  printed  values. 


The  format  is  the  same  as  for  *IV 


is  either  a  constant  alone  or  Pi  =  constant.  Values  are 
separated  by  blanks.  If  a  constant  appears  alone,  it  is 
assumed  to  be  the  value  of  the  next  PAR(i).  If  no  value 
is  given  for  a  PAR(i),  it  is  assumed  to  be  zero. 

$  is  a  required  terminator 


EXAMPLES: 

1.  *PAR  P2=6.  75  4.0  P5=8.2  $ 

This  gives:  PAR(1)=0,0 

PAR(2)=6.  75 
PAR  (3)  =4.  0 
PAR(4)=0.  0 
PAR(5)=8.  2 

2.  Use  of  the  PAR  array  to  change  parameters  easily  in  a  series  of  pro¬ 
gram  runs  is  illustrated  in  the  following  example: 

If  one  of  the  derivative  equations  is 

DY(3)=SIN(15.)  *  Y (1 )  +  C0S(2O.)*Y(2)  +  SIN(25.  )*Y (3) 

and  it  is  desired  to  vary  the  angles  over  a  range  of  values,  the 
derivative  equation  can  be  written 

DY(3)=SIN(PAR(1))*Y(1)  +  C0S  (PAR  (2)  )*  Y  (2)  +  SIN(PAR  (3  ))*Y(3) 

and  the  following  *PAR  cards  used,  one  per  run,  without  changing 
the  equation 


PAR 

15.  0 

20.  0 

25.  0 

$ 

PAR 

30.  0 

35.  0 

40.0 

$ 

PAR 

72.  0 

75.  0 

78.0 

$ 

where 

par. 


’I 


*PAR  and  *IV  cards  are  particularly  useful  in  making  multiple  runs, 
that  is,  in  making  several  runs  of  the  same  set  of  equations  as  One  job, 
varying  only  the  initial  conditions,  PAR's,  or  perhaps  run  times  from  one 
run  to  the  next  (see  Appendix  F-l,  "Multiple  Runs").  Note  that  since  PAR's 
and  IV's  retain  their  values  from  run  to  run  within  a  job  until  they  are  reset 
by  the  user,  only  those  being  changed  for  a  particular  run  need  to  be  redefined. 

D.  INITIAL  CONDITIONS  AND  INPUTS  TO  BE 

COMPUTED  (*ICC0MP. .  .  *ENDIC) 

If  any  computation  is  necessary  to  set  initial  conditions  or  parameter 
values,  the  user  will  need  to  write  an  initial  computations  section  into  his 
coding.  The  start  of  this  section  is  signaled  by  the  card  *ICC0MP,  and  the 
end  is  signaled  by  *ENDIC.  Between  these  two  cards,  standard  F0RTRAN 
and  WHELP  statements  may  be  used,  and  their  sequence  is  governed  only 
by  the  usual  rules  of  FORTRAN  and  WHELP.  All  statements  in  this  section 
will  become  part  of  the  initial  computations  subroutine,  ICC0MP,  written  by 
the  ESP  precompiler  program,  PREC0MP. 

This  subroutine  is  executed  once  (and  only  once)  per  run,  so  it  is  the 
proper  place  to  perform  any  operation  that  is  to  be  done  only  once  before 
starting  the  solution,  such  as  computing  PAR's,  computing  initial  values, 
reading  in  data,  or  rewinding  tapes. 

*ICC0MP.  .  .  .  *ENDIC  may  be  used  in  addition  to  or  in  place  of  *PAR 
and  *IV  cards.  Since  the  resulting  subroutine  is  executed  after  *IV,  *PAR, 
and  *RUN  have  been  encountered  and  processed,  the  user  may  input  constants 
via  *PAR  or  -TV  and  then  safely  use  them  on  the  right-hand  side  of  expres¬ 
sions  within  *ICC0MP, 

The  program  segment  *ICC(t)MP.  .  .  *ENDIC  should  be  positioned  after 
all  user  supplied  routines,  but  before  all  run  time  cards  such  as  *IV,  *PAR, 
or  *RUN  (refer  to  Appendix  A-3  for  deck  structure). 
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EXAMPLE: 


*ICC0MP 

DIMENSION  R  (6) 

EQUIVA LENCE(PAR(1),  R),  (PAR(7),  RD),  (PAR(8),  DR) 
DATA  PI/3.  14159/ 

RD=180.  /PI 
DR  =  PI/  1 80. 

C  INPUT  R  ANGLES  AS  PARS  IN  DEGREES  AND  CONVERT  HERE. 
DO  5  1=1,6 
R(I)  =  DR*PAR(I) 

5  CONTINUE 
*ENDIC 

*  PA  R  10.  20.  30.  40.  50.  60.  $ 

*RUN  3  0.  0.5  10.0  $ 


DATA  INPUT  FROM  CARDS  OR  USER  FILES 


Because  an  ESP  job  sometimes  requires  the  input  of  data  which  cannot 
be  conveniently  handled  by  *PAR  cards  or  DATA  statements  in  ICC0MP,  the 
user  also  has  the  option  of  reading  it  in  with  READ  or  NAMELIST  statements. 
The  data  itself  may  either  be  on  cards  as  part  of  the  user's  deck  or  it  may 
be  in  the  form  of  a  user-created  file  (stored  on  magnetic  tape  or  disc 
storage).  In  either  case,  the  best  place  to  read  it  from  is  within  *ICC0MP.  .  . 
*ENDIC. 

If  the  data  exists  on  a  manageable  number  of  cards,  the  data  cards  may 
be  placed  immediately  after  the  *RUN  card(s)  for  the  case(s)  to  which  they 
apply.  Since  PREC0MP  copies  all  cards  beginning  with  the  first  run  time 
control  card  onto  TAPE  12,  data  cards  so  introduced  will  exist  there  during 
job  execution  along  with  the  normal  run  time  control  cards.  These  data 
cards  may  be  read  as  if  they  were  on  the  standard  INPUT  file,  e.g.,  READ  100, 
A,  B,  C. 

NOTE 

The  user  must  take  care  to  read  exactly  the  correct 
number  of  cards  since  proper  execution  of  *GRAPH  or 
any  other  run  time  cards  depends  on  the  proper  position¬ 
ing  of  the  input  file.  Also,  input  data  may  occupy  only 
the  first  72  columns,  as  columns  73-80  will  not  be  copied 
to  TAPE  12. 


\.y * r  " 


Sometimes,  however,  it  may  be  more  convenient  to  read  input  data 
from  the  user's  own  file.  This  would  be  true,  for  example,  if  the  amount  of 
data  is  very  large,  if  it  is  necessary  to  test  for  end  of  file  to  terminate 
reading  data  for  a  case,  or  if  t! —  user  already  has  the  data  on  a  stored  file 
of  some  sort.  If  this  is  the  case,  the  user  should: 

•  Create  and  save  the  data  file  (if  it  does  not  already  exist),  being 
sure  to  write  in  End-of-Files  as  he  will  need  them. 

•  Add  the  proper  job  control  cards  to  assign  the  saved  file  to  the 
job  and  give  it  a  logical  name,  say  TAPEn. 

•  Write  his  own  MAIN  program  (refer  to  Appendix  D-3-b),  adding 
TAPEn  to  the  files  declared  on  the  PROGRAM  card. 

•  Read  the  data  from  the  file  in  ICCOMP  by  using 

READ  (n,  format)  list  or  [  formatted  read] 

READ  (n,  name)  [namelist  read] 

READ  (n)  list  [unformatted  read] 

For  more  on  data  input  used  with  stacked  or  multiple  runs,  refer 
to  Appendix  F. 

F.  INPUTS  TO  CONTROL  ACCURACY 

In  addition  to  the  inputs  described  above,  there  are  several  optional 
inputs  which  may  be  used  to  control  the  accuracy  of  the  solution  and  the  timing 
accuracy  of  discontinuities.  All  have  default  values  but  may  be  user-defined 
by  means  of  special  control  cards  placed  in  any  order,  among  the  run-time 
data  cards,  that  is,  after  the  ^DERIVS,  ^OUTPUT,  and  *ICCC5MP  sections, 
but  before  *RUN.  *EPS  and  *Q  control  the  solution  accuracy  directly  and  are 
described  in  Section  IV-A-4.  #HSW,  *HSWM,  and  *HSWE  are  used  to  control 
the  allowable  timing  error  in  SWTCH's,  SWMEM's,  and  EVENT'S,  respectively, 
and  are  discussed  in  Section  V-F. 

G.  MISCELLANEOUS  INPUTS 

Several  other  special  input  cards  will  be  read  and  interpreted  by 
PREC0MP.  Their  use  is  optional  and  they  are  provided  mainly  for  the  con¬ 
venience  of  the  user. 

•  t  •  t  f 
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1 .  Print  Headings 

The  print  labels  which  appear  on  automatically  formatted  output  are 
normally  specified  on  the  *PRINT  card  (refer  to  Section  VI-A)  from  which 
they  are  read  by  PREC0MP  and  written  onto  a  run-time  data  card  called 
-^HEADINGS,  which  the  user  will  notice  is  printed  with  the  other  run  time 
data  cards  at  the  beginning  of  his  output.  The  user  may,  however,  supply 
the  ^HEADINGS  card  himself,  in  which  case  his  card  completely  supersedes 
the  card  written  by  PREC0MP.  The  number  of  labels  on  this  card  is  the 
number  of  output  variables  which  will  be  printed,  even  if  it  is  less  than  the 
number  of  variables  specified  on  the  print  card.  (Thus  a  ^HEADINGS  card 
with  no  labels  can  be  used  to  suppress  printout.)  The  format  is 


where 

SHEADINGS  starts  in  column  1 

label.  10  hollerith  characters  (8  on  IBM)  with  no  embedded 

blanks,  which  will  be  used  to  label  the  ith  output  variable. 

n  is  the  number  of  variables  that  will  be  printed.  It  should 

be  the  same  or  less  than  the  number  of  variables  listed  on 
the  *PRINT  card. 

$  is  a  required  terminator 

2.  Title  for  Printed  Output 

Placing  the  *TITLE  card  among  the  run-time  data  cards,  somewhere 
before  the  *RUN  card,  causes  whatever  is  in  columns  8-71  to  be  used  as  the 
title  printed  at  the  top  of  every  page  of  output.  The  format  is 

f - 1 


3. 


Program  Control 

The  following  two  cards  may  be  used  anywhere  among  the  run-time  data 
cards,  to  facilitate  program  control: 


-RETURN 


-STOP 


This  causes  program  control  to  return  to 
PROGRAM  MAIN  at  the  point  it  is  en¬ 
countered,  typically  following  a  -RUN  or 
-GRAPH  command.  It  thereby  permits 
execution  of  other  statements  in  MAIN, 
such  as  calls  to  other  subroutines  or 
printing  which  is  to  be  done  only  at  the  end 
of  a  run  (see  Section  VI-C-5). 


This  causes  job  termination  at  the  point  it  is 
encountered.  It  should  be  used  after  *RUN  if 
no  -GRAPH  follows  and  normal  termination  is 
desired  at  that  point. 
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DECK  STRUCTURE 


APPENDIX  A 


CARD  FORMATS  AND  DECK  STRUCTURE 


A  -  1.  GENERAL  FORMAT  RULES 


General  rules  regarding  ESP  card  formats  are  the  following: 

•  All  ESP  "shorthand"  cards  begin  with  an  asterisk  (*)  in 
column  1  and  the  first  letter  of  the  key  word  in  column  2. 

•  Items  on  cards  are  separated  either  by  blanks  or  by  $ 
depending  on  their  nature.  (See  specific  cards.) 

•  Fixed  length  cards  generally  require  no  terminator,  but 
those  of  variable  length  are  terminated  by  a  required  $. 

•  FORTRAN  statements  begin  in  column  7  and  follow  the  usual 
rules  of  F0RTRAN. 

•  WHELP  statements  follow  the  usual  format  of  WHELP 
(Appendix  H). 

General  rules  regarding  use  of  comment  cards: 

•  Comment  cards  are  denoted  only  by  a  C  in  column  1. 

They  may  appear  anywhere  within: 

•  user  (sub)programs 

•  the  derivative  computation  section,  delimited  by  *DERIVS 
and  *ENDDERIVS 

•  the  output  routine  delimited  by  *0UTPUT  and  -'END0UT 

•  the  initial  computation  routine  delimited  by  -TCC0MP 
and  *ENDIC 

•  the  run-time  data  cards  (provided  they  are  between 
command  cards). 

•  Comment  cards  may  not  be  used: 

•  within  the  data  picked  up  in  response  to  a  command  card, 
e.g.,  within  the  *PRINT  specification. 

•  between  any  routines,  i.e.,  user  routines,  after 
*ENDDERIVS,  *END0UT  or  *ENDIC,  or  after  -’'-PRINT 

^  whet^  not  within  0UTPUT. 
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A -2.  SUMMARY  OF  CARD  FORMATS 


Section 

Reference 


•  Defining  derivatives  and  discontinuities: 

*DERIVS 

*ENDDERIVS 

*B  LOCK  1  0  Y (i)  e.n  $ 

*BLOCK  2  ax  aQ  01  0Q  Y (i)  Y(j)  e.n  $ 

*SWTCH  i  0+  $  0  $  control.  $ 

*SWMEM  i  inpuL  $ 

•  Defining  output  and  initial  computations: 

*<2UTPUT 

*END0UT 

*PRINT  labelj=expressionj  $.  .  .  label^expression^  $  $ 

*ICC0MP 

*E  NDIC 

•  Run-time  data  inputs: 

*IV  iv.  iv_ . . . iv  $ 
it.  n 

-:-PAR  px  P2-..Pn  $ 

*S  WITCHES  n 
*SWMEMCNT  n 
*NEVENT  n 

*SWMEMSET  n:  . .  n£ 

*SWMEMDATA 


C2  •  *  •  c 10 


1  c . 
n  1 


*•  c10  $ 


III-C 

III-C 


VI- A  -  2 


VI- A  -  1 


VII- D 


VII-  B 
VII- C 


V-D-l-b 


V-B-4 

V-B-3 


Section 

Reference 


♦  MAX PLOTS  n 

♦RUN  neq  hprint^  tf^  hprint^  tf^  .  .  .  $ 
•  Accuracy  control: 

*EPS  e 


♦Q 

ql  q2 

qn 

$ 

♦  HSW 

hl  h2 

h 

n 

$ 

♦  HSWM 

hl  h2 

h 

n 

$ 

♦HSWE 

hl  h2  - * 

h 

n 

$ 

•  Miscellaneous  inputs: 

♦  TITLE - title - 

♦HEADINGS  label,  label.,  ...  Label  $ 
12  n 

♦STOP 

♦RETURN 

RK2 

♦METHOD  RK4 
PC 

•  Plotting: 

♦GRAPH  n  n  [size][grid][ scaling][type] 

x  y 

[title] 

[X  title] 

[Y  title] 


VI-  B 

VII-  A 

IV-A-4 

IV- A-4 

V- F 
V-F 
V-F 

VII-  G-  2 
VII-G-1 
VII-G-3 
VII-G-3 

IV- A-  1 


VI-B-2 


A  -  3 


A  -  3.  USER'S  DECK  STRUCTURE 


! 

% 

f 

s 

t 


(CDC) 

[job  control  cards 
7-8-9  card 


[*METH0D 


RK2 

RK4 

PC 


(IBM) 

J  C  L  cards 
//SYSIN  DD  * 


[PR0GRAM  MAIN  (if  written  by  user)  ^PROGRAM  ON  IBM 

All  user-supplied  subroutines,  including  ICC0MP,  0UTPUT, 
SWINPT,  SWMEMN,  EVENTS,  and  N0TIFY,  if  the  user  pro¬ 
vides  the  entire  routine.  (See  Appendix  C-4  reserved  names.) 

Coding  which  defines  derivatives  and  discontinuities,  in¬ 
cluding  all  *SWTCH,  *SWMEM,  and  -:BL0CK  statements 
used,  beginning  with  *DERIVS  and  ending  with  *ENDDERIVS. 

[Coding  which  defines  the  output:  either  *PRINT  or 
*0UTPUT.  .  .  *END0UT. 

[Coding  which  defines  the  initial  computations :  -ICC0MP.  .  . 
*ENDIC. 

Run-time  data  cards,  in  any  order,  including  *IV,  *PAR, 
-SWITCHES,  -SWMEMCNT,  -SWMEMDAT A,  -SWMEMSET, 
.-EPS,  -Q,  -HSW,  -HSWM,  -HSWE,  *NEVENT,  *MAXPL0TS 

[  -  RUN 

[Any  input  data  to  be  read  by  a  READ  format,  list  or 
READ  namelist. 


[  -GRAPH 


[  -ST0P  or  -RETURN 


6-7-8-9 


/- 

JC L  cards 

/- 


Only  those  sections  or  cards  marked  by  an  -  are  required.  All  others  are 
optional. 

^The  order  of  these  three  segments,  *DERIVS,  -OUTPUT,  and  *ICC0MP, 
is  interchangeable. 
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CONTROL  CARDS  AND  FILE  USAGE 


B-l.  CDC  CONTROL  CARDS 


(SCOPE  2.  15  OPERATING  SYSTEM) 


NOTE 

The  following  control  card  examples  apply  to  the 
operating  system  in  use  at  the  publication  date  of 
this  manual.  Subsequent  changes  in  operating  sys¬ 
tem  or  control  cards  required  will  be  documented  as 
they  occur. 


ESP  USED  WITH  WHELP 


$PGMR.  . .  . 

$PARAM.  . .  . 

ATTACH,  LIB  I,  2NEWRESP. 
ATTACH,  LIB 2,  3FTNPL0TLIB. 
LIBRARY,  LIB  1,  LIB2. 
PRECOMP. 

WHELP,  IMSORC. 

FTN,  I=TA PE15. 

LGO. 

ESP  USED  WITHOUT  WHELP 

$PGMR.  . . . 

$PARAM.  . . . 

ATTACH,  LIB  I,  2NEWRESP. 
ATTACH,  LIB2,  3  FTN  PLOT  LIB. 
LIBRARY,  LIB1,  LIB2. 
PRECOMP. 

FTN,  I=IMSCRC. 

LG  0. 


To  get  hard  copy  of  film  plots,  add  after  LGO.  :  HARDCPY. 
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B-2.  IBM  CONTROL  CARDS  (JCL) 


ESP 


00010 

00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 


00120 


00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 

00220 

00230 

00240 

00250 

00260 

00270 


(IBM  3033  MVS  OPERATING  SYSTEM) 
WITH  WHELP 


/  /  Z 1 85  JOB... 

//  MSG  LEVEL.  .. 

/*J0BPARM  ACCT... 

//  EXEC  PGM=PREC0MP 

//STEPLIB  DD  DSNr#4606.  ESP.  LIB  (PREC0MP),  DISP=SHR 
/  /SYSPR1NT  DD  SYS0UT=A 

/  / T A  PE  1 1  DD  DSN  =  &TA  PE  1 1,  DISP=(NEW,  PASS),  UNIT  =  VI0, 

//  SPACE  =  (TRK,  (10.  5)).  DCB  =  (RECFM  =  FB.  LRECL=80.  BLKS1ZE=800) 

/  /  T  A  PE  1 2  DD  DSN  =  &TAPE12,  DISP=(NEW,  PASS).  UNIT=VIO, 

//  SPACE  =  (TRK,  (10,  5)),  DCB  =  (RECFM  =  FB.  LRECL=80,  BLKSIZE=800) 

/  /  TAPE  14  DD  DSISI=&<&.T^PE14.  DISP=(NEW,  DELETE),  UNIT  =  VlO, 

//  SPACE  =  (TRK,  (10*  5)),  DCB  =  (RECFM  =  FB,  LRECL=80,  BLKSIZE  =  800) 

//TAPE16  DD  DSN  =  &,&TAPE16.  DISP=(NEW,  DELETE).  UNIT  =  VI0, 

//  SPACE  =  (TRK,  (10.  5)),  DCB  =  (RECFM=FB,  LRECL=80,  BLKS1ZE  =  800) 

//SYSIN  DD  DSN=#USERID.  FILENAME.  DATA.  DISP-SHR 
(if  user  input  resides  on  a  permanent  file) 
or 

//SYSIN  DD  * 

(ESP  source  program  cards) 

//  EXEC  PGM = WHELP 

//STEPLIB  DD  DSN=#4606.  ESP.  LIB  (WHELP),  DISP=SHR 
//SYSPRINT  DD  SYS0UT=A 

//SYSIN  DD  DSN=&TAPE11,  DISP=(©LD,  DELETE) 

/  /  TAPE  15  DD  DSN=8tTAPE15,  DISP*(NEW,  PASS),  UNIT  =  VI0, 

//  SPACE=(TRK.  (10,  5)).  DCB  =  (RECFM  =  FB,  LRECL=80,  BLKSIZE=800) 

/ /  TAPE  1 1  DD  DSN=kTEMP,  DISP=(NEW,  DELETE) ,  UNIT=VI0, 

//  SPACE  =  (TRK,  (10,  4)),  DCB  =  (RECFM  =  FB,  LRECL=80,  BLKSIZE=800) 

//  EXEC  F0RTXCLG,  CPARM  =  'N0F0RMAT,  AD(DBL),  MAP',  C0ND.  LKED=EVEN, 
//  C0ND.GO=EVEN,  LPARM  =  LET 

//FORT. SYSIN  DD  DSN=&TAPE15,  DISP=(0LD,  DELETE) 

/  /LKED.  SYS  LIB  DD  DSN=0PUS.  P077.SUBLIB,  DISP=SHR 

//GO.  FT11F001  DD  DSN  =  fc& TAPE11 ,  DISP=(NEW,  DELETE),  UNIT=VI0, 

//  SPACE  =  (TRK,  (10,  5)),  DCB  =  (RECFM=VBS,  B LKSIZE=6440,  LRECL=16004) 

//GO.  FT12F001  DD  DSN  =  &TAPE12,  DISP=(0LD,  DELETE) 

/* 
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ESP 


00010 

00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 


00120 


00210 

00220 

00230 

00240 

00250 

00260 

00270 


WITHOUT  WHELP 


'  Z 1 85  JOB... 

;  ,  MSG  LEVEL. .. 

/*J0BPARM  ACCT... 

//  EXEC  PGM  =  PREC0MP 

//STEPLIB  DD  DSN=#4606.  ESP.  LIB(PRECOMP),  DISP=SHR 
//SYSPRINT  DD  SYS0UT=A 

//TAPE11  DD  DSN=&TAPE11,  DISP=(NEW,  PASS),  UNIT  =  V10, 

//  SPACE=  (TRK,  (10,  5)).  DCB  =  (REC  FM=  FB,  LRECL=80,  BLKS1ZE=800) 

/  /  T  A  PE  12  DD  DSN=&TAPE12,  DISP=(NEW.  PASS).  UNIT  =  V10, 

//  SPACE  =  (TRK,  (10,  5)),  DCB  =  (REC FM=FB,  LRECL=80,  BXJ<SIZE=800) 

/ /  TAPE  14  DD  DSN=& &TAPEI4,  DISP  =  (NEW,  DELETE),  UNIT  =  VI0, 

//  SPACEs (TRK,  (10,  5)),  DCB=(RECFM  =  FB,  LREC L=80) ,  B LKS1ZE=800) 

/ /TAPE  16  DD  DSN=&&iTAPE16,  DISP=(NEW,  DELETE),  UNIT=VI0, 

//  SPACE= (TRK,  (10,  5)),  DCB=(RECFM  =  FB,  LRECL=80,  BLKSIZE=800) 

//SYSIN  DD  DSN=#USERID.  FILENAME.  DATA,  DISPsSHR 
(if  user  input  resides  on  a  permanent  file) 
or 

//SYSIN  DD  * 

(ESP  source  program  cards) 

//  TEXEC  F0RTXCLG,  CPARM='NOFORMAT,  AD(DBL),MAP',C0ND.  LKED=£VEn 
//  C0ND.G0  =  EVEN,  LPARM  =  LET 

//F0RT. SYSIN  DD  DSN=&TAPE11,  DISP=(0LD,  DELETE) 

//LKED.SYSLIB  DD  DSNs0PUS.  P077.  SUB  LIB,  DISP=SHR 

//GO.  FT11F001  DD  DSN=&&TAPE11,  DISP=(NEW,  DELETE),  UNIT=VI0, 

II  SPACE=(TRK,  (10,5)),DCB  =  (RECFM=VBS,BLKS1ZE*6440,  LRECL=lfe004) 

/ / G0.  FT12F001  DD  DSN=&TAPE12,  DISP=(0LD,  DELETE) 

/* 
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B -  3.  FILE  USAGE  FOR  ESP  WITHOUT  WHELP 


B-4.  FILE  USAGE  FOR  ESP  WITH  WHELP 
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PROGRAM  VARIABLES  AND  RESERVED  NAMES 


C-l.  VARIABLES  PASSED  THROUGH  CA  LLTNG  SEQUENCES 


Variable 

Routines  to  which 
it  is  passed: 

DY  (100) 

The  derivative  array 

DERIVS 

OUTPUT 

IE  VENT 

An  integer  indicating  the  number  of  the 
EVENT  being  reported  to  SUBROUTINE 
NOTIFY 

NOTIFY 

PL0T(1OO) 

An  array  for  storing  the  current,  va lue  (s ) 
of  the  plotted  variables  (Equivalent  to 
VPLOT  in  COMMON  BL0CK  UNIP1) 

OUTPUT 

PRINT  (60) 

An  array  for  storing  the  current  value(s) 
of  the  printed  variables 

OUTPUT 

ST0P 

A  variable  which  stops  the  current  run  if 
nonzero  (Equivalent  to  T  >  TFINAL) 

DERIVS 

OUTPUT 

T 

The  independent  variable,  usually  time 

DERIVS 

OUTPUT 

ICC0MP 

SWMEMN 

SWINPT 

VA  LUES(50) 

An  array  for  storing  the  inputs  to  SWTCHs, 
SWMEMs,  or  EVENTs 

SWINPT 

SWMEMN 

EVENTS 

Y  (100) 

The  dependent  variable  array 

DERIVS 

OUTPUT 

iccOmp 

SWINPT 

SWMEMN 

C-l 


C-2.  VARIABLES  PASSED  THROUGH  COMMON  BLOCKS 


Variable 

Block  Name 

Variable 

BUFFER  (80) 

READIN 

NDISPR 

C0LCNT 

READIN 

NEQ 

C0NSTS  (50,  10) 

SWHPAR 

NEVENT 

DY  (100) 

BASIC 

NFIRST 

DY (100,  9) 

BLANK 

NHEAD 

EPS 

MISCEL 

NL0CAL 

FIRSTP 

RKC0NT 

N0PL0T 

FIXSTP 

STPC0N 

NPAGE 

H 

STPC0N 

NP0INT 

HEAD(60) 

UNIP2 

NTAP11 

HHMAX 

HMAXMN 

NUMSTP 

HHMIN 

HMAXMN 

0UT  (60) 

HMAX 

STPC0N 

PARflOO) 

HMIN 

STPC0N 

P  L0T  (2000) 

HP 

STPC0N 

Q(100) 

HSW  (50) 

MISCEL 

ST0P(1) 

HSWE  (50) 

MISCEL 

SCRl(200) 

HSWM(50) 

MISCEL 

SWDBUG 

IF0RM(3) 

UNIP2 

SWMEM(50,  4) 

ISWTYP 

SWHPAR 

SWSET  (50) 

J  LINE 

UNIP2 

SWTCH(50) 

JSTART 

STFPAR 

TO 

KFLAG 

STFPAR 

TERMCH 

KSV 

SWHPAR 

TF 

LINES 

UNIP2 

TIT  LE  (8) 

MF 

STFPAR 

T0DA  Y 

MAX 

UNIPI 

TP 

MAXCHR 

READIN 

VA  LE  VS(50,  2) 

MAXC0L 

READIN 

VA  LMEM(50,  2) 
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Block  Name 

NDISPR 

MISCEL 

SWT CHS 

READIN 

UNIP2 

UNIP1 

UNIP1 

UNIP2 

UNIPI 

HMAXMN 

HMAXMN 

UNIP2 

PARS 

BLANK 

MISCEL 

MISCEL 

BLANK 

SWDBUG 

SWTCHS 

SWHPAR 

SWTCHS 

BASIC 

READIN 

BASIC 

UNIP2 

UNIP2 

BASIC 

MISCEL 

MISCEL 


Variable 


Block  Name 


Variable 


Block  Name 


MAXDER 

STFPAR 

VA  LUES  (50,  2) 

MISCEL 

MAXMEM 

SWTCHS 

VPL0T(1OO) 

UNIP1 

MAXSWS 

SWTCHS 

Y(100,  9) 

BLANK 

MX  L 

UNIP2 

Y0(100) 

BASIC 

NALTER 

SWHPAR 

YPRNT(IOO) 

BASIC 

NCHNG 

SWHPAR 

C-3.  ALPHABETICAL  LIST 

OF  COMMON  BLOCKS  AND 

THEIR 

CONTENTS 

NOTE 

Blocks  SWT CHS  and  PARS  are  written  by  PREC0MP  as 
part  of  SUBROUTINES  DERIVS,  OUTPUT,  ICC0MP, 
SWINPT,  and  SWMEMN.  All  other  blocks  must  be  included 
by  the  user  if  he  wishes  to  reference  them.  Block  lengths 
are  given  in  octal  words  for  CDC  use  and  hexadecimal  bytes 
for  IBM  use  since  these  are  the  bases  used  to  list  block  size 
on  maps  generated  by  the  two  computers.  Decimal  size  of 
each  block  can  be  easily  obtained  from  the  dimensions  given 
with  each  variable  name. 


■ 

Block 

Length 

CDC 

(Octal 

IBM 

(Hex) 

C0MM0N/ BASIC  /  TO,  TF,  TP.  Y0(100).  YPRNT(IOO),  DY(100) 

457 

978 

TO 

The  initial  value  of  the  independent 
variable 

1 

8 

TF 

The  current  final  value  of  the  independent 
variable  (May  be  changed  by  the  user's 
program) 

1 

8 

TP 

The  last  print  time 

1 

8 
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Block  Length 


CDC 

(Octal) 

IBM 

(Hex) 

Y0(100) 

The  initial  conditions  to  be  used  when 
(next)  *RUN  is  encountered.  Set  =  0 
in  ESP. 

144 

320 

YPRNT  (100) 

The  value(s)  of  the  independent  variables 
at  the  last  print  time.  After  a  run  this 
contains  the  "final"  values  of  the  Y's. 

144 

320 

DY(100) 

The  value(s)  of  the  derivatives  at  the  last 
print  time.  After  a  run  this  contains  the 
"final"  values  of  the  DY's. 

144 

320 

C0MM0N/BLANK/PL0T  (2000),  SCR1  (200) ,  Y(100,  9) . 

DY (100,  9)f 

7640 

7D00 

PLOT  (2000) 

Plot  buffer 

3720 

3E80 

SCR1  (200) 

Used  internally  as  working  space 

310 

640 

Y (100,  9) 

Y  array,  including  past  values  of  Y's 

1604 

1C20 

DY(100,  9) 

DY  array,  including  past  values  of  DY1  s 

1604 

1C20 

C0MM0N/HMAXMN/HHMAX.HHMIN,  NUMSTP,  NTAP11 

4 

18 

HHMAX 

The  maximum  stepsize  used  thus  far  in 
the  run.  Automatically  printed  at  end  of 
run.  May  be  tested  or  printed  by  user, 
but  not  changed. 

1 

8 

HHMIN 

The  minimum  stepsize  used  thus  far  in 
the  run.  Automatically  printed  at  end  of 
run.  May  be  tested  or  printed  by  user 
but  not  changed. 

1 

8 

NUMSTP 

The  number  of  integration  steps  taken  thus 
far  in  the  run.  Automatically  printed  at 
end  of  run.  May  be  tested  or  printed  by 
user  but  not  changed. 

1 

4 

NTAP11 

The  number  of  data  frames  written  onto 
TAPE11  for  plotting.  It  is  automatically 
printed  at  the  end  of  the  run. 

1 

4 

The  storage  in  this  common  block  is  used  differently  by  different  subroutines, 
but  this  describes  the  most  common  and  generally  relevant  use.  The  user 
may  wish  at  times  to  know  the  contents  of  this  block,  but  should  not  alter  them. 
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C0MM0N/MISCEL/ST0P(1),  Q(iOO),  EPS,  HSW(50),  HSWM(50), 
HSWE(50),  VALUES(50,  2),  VA  LMEM(50,  2),  VALEVS(50,  2). 
NEQ 


ST0P(1) 


Q  (100) 


HSW  (50) 


HSWM  (50) 


HSWE  (50) 


VA  LUES  (50,  2) 


A  variable  which  stops  the  current  run  if 
nonzero  (equivalent  to  T  >  TFINAL) 

Q(i)  is  used  to  compute  a  maximum  allow¬ 
able  absolute  error  in  Y(i).  It  is  set 
dynamically  to  MAX(Q(i),  |  Y(i)J  ). 

EPS  is  used  to  compute  relative  error  in 
Y  (i) : 


error  in  Y(i) 

Q(i) 


HSW(i)  is  the  allowable  timing  error  in 
determining  SWTCH(i),  normally  set  on 
the  *HSW  data  card. 

HSWM(i)  is  the  allowable  timing  error  in 
determining  SWMEM(i),  normally  set  on 
the  *HSWM  data  card. 

HSWE(i)  is  the  allowable  timing  error  in 
determining  EVENT  (i),  normally  set  on 
the  *HSWE  data  card. 

VALUES(i,  1)  and  VALUES(i,  2)  store  the 
current  and  previous  values  of  the  inputs 
to  *SWTCH(i),  alternately. 


VALMEM(50,  2)  VALMEM(i,  1)  and  VA LMEM(i,  2)  store  the 
current  and  previous  values  of  the  inputs 
to  *SWMEM(i),  alternately. 

VALE  VS  (50,  2)  VALEVS(i,  1)  and  VALEVS(i,  2)  store  the 

current  and  previous  values  determining 
EVENT (i),  alternately. 


NEQ 


The  number  of  derivative  equations  to  be 
integrated:  set  by  the  user  on  the  *RUN 
card. 


Block  Length 


CDC  IBM 
(Octal)  (Hex) 


1051  1144 


144  320 
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Block  Length 

CDC  IBM 
(Octal)  (Hex) 


C0MM0N  /  NDISPR  /NDISPR  1 

NDISPR  A  flag  which  controls  printing  at  switching  1 

points.  If  0,  no  print  at  switch  times;  if  1, 


one  print  at  switch  times;  if  2,  print  and 
plotting  occurs  on  left  and  right  of  each 
switch.  (NDISPR  =  1,  nominally) 

C0MMpN/PARS/PAR(lOO)  144 


PAR(IOO)  An  array  for  storing  and  automatically  144 

transmitting  user  variables,  which  may  be 
easily  input  on  a  *PAR  card. 

C0MM0N/ REA  DIN/ C0LCNT  ,  BUFFER  (80),  MAXC0L, 

MAXCHR,  NFIRST,  TERMCH  12  5 


C0L.CNT  Used  internally  by  READIT:  points  tc  1 

beginning  of  next  scan. 

BUFFER(80)  Used  internally  by  READIT:  the  current  120 

card  in  80A1  FORMAT. 

MAXC0L,  Used  internally  by  READIT:  the  last  1 

column  to  be  scanned. 

MAXCHR  Used  internally  by  READIT:  maximum  1 

number  of  characters  to  be  picked  up 
in  a  hollerith  field 

NFIRST  Used  internally  by  READIT:  points  to  1 

beginning  of  field  just  read. 

TERMCH  Used  internally  by  READIT:  any  hollerith  1 

character  to  mark  the  end  of  a  field. 

C0MM0N/  RKC0NT  /  FIRSTP  1 


FIRSTP  A  flag  set  to  1.0  by  routine  ESPCTLif  1 


Runge-Kutta  is  used  to  indicate  the  beginning 
of  each  step  in  fixed  step  mode  or  the 
beginning  of  each  pair  of  steps  in  the  variable 
step  mode.  Otherwise,  FIRSTP=0. 


4 

4 


320 

320 


158 

4 

140 

4 

4 


4 

8 

8 

8 
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Block  Length 


C0MM0N/STFPAR/MF,  KF BAG,  JSTART,  MAXDER 


CDC  IBM 

(Octal)  (Hex) 

4  10 


MF 

KFLAG 


JSTART 


MAXDER 


Used  internally 

A  flag  returned  from  the  integration 
routines  to  indicate  success  or  failure 
of  the  integration  step  just  taken. 

KFLAG  =  1  indicates  error  exceeded 
bounds  and  a  warning  message  willbe 
printed. 

A  flag  used  to  indicate  the  start  (restart) 
or  continuation  of  integration.  JSTART  =  0 
when  integration  is  starting  or  restarting. 
JSTART  =  1  when  integration  is  continuing 
on  from  previous  steps. 

Used  internally 


C0MM0N/STPC0N/HP.H,  FIXSTP,  HMIN,  HMAX 


HP 


H 


FIXSTP 


HMIN 

HMAX 


The  current  printing  interval.  This  is 
normally  changed  from  the  *RUN  card 
but  may  be  changed  by  the  user's  pro¬ 
gram  during  the  run  and  must  be  >  0. 

The  current  integration  stepsize.  If 
set  ^  0  in  ICC0MP  this  H  will  be  tried 
first. 

The  actual  stepsize  selected  by  the  user 
for  fixed  stepsize  integration  using  all 
Runge  -  Kutta 

A  lower  limit  on  the  stepsize,  nominally  0. 
May  be  set  by  user. 

The  maximum  stepsize  permitted, 
nominally  1. 0E50.  May  be  set  by  user. 


1 

1 


1 


1 

5 

1 


1 


1 


1 

1 


4 

4 


4 


4 

28 

8 


8 


8 


8 

8 
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Block  Length 


CDC 

(Octal) 

C0MMQN/SWDBUG/SWDBUG  1 

SWDBUG  Logical  variable  which  controls  printing  of  1 

data  for  switch  debugging.  If  SWDBUG  = 

.FALSE,  (default)  no  print.  If  SWDBUG  = 

.TRUE,  print  data  to  aid  in  debugging  of 
switches . 


C0MM0N/SWHPAR/NCHNG,  NALTER,  ISWTYP,  KSV, 

C0NSTS(5O,  10),  SWSET (50)  1052 


NCHNG 


NALTER 


A  flag  indicating  whether  any  switches  1 

have  just  changed  state  during  an  integra¬ 
tion  step 

A  flag  used  internally  by  SWTCHE  1 


ISWTYP 


Used  internally 


1 


KSV 


Used  internally 


1 


C0NSTS(5O,  10)  C0NSTS(i,  j)  is  the  constant  Cj  for  764 

SWMEMi.  Although  C0NSTS  are  normally 
defined  on  the  *SWMEMDATA  card,  the 
user  may  include  common  block  SWHPAR 
and  define  the  CONSTS  in  ICC0MP.  (No 
error  test  is  made  on  C0NSTS  so  defined.  ) 


SWSET (50)  The  vector  of  values  used  to  initialize  62 

SWMEMS  in  saturation  rather  than  dead¬ 
band.  Normally  input  on  the  *SWMEMSET 
card. 


C0MM0N/SWTCHS/SWTCH(5Q),  SWMEM(50,  4),  MAXSWS, 

MAXMEM,  NEVE  NT  375 

SWTCH(50)  The  magnitude  is  1+  the  number  of  times  62 

SWTCHi  has  switched.  The  sign  is  the 
current  sign  of  the  input.  On  the  first  call 
to  DERIVS  following  a  switching,  all  switches 
which  have  changed  state  have  their  magni¬ 
tudes  increased  by  0. 5. 


IBM 

(Hex) 

4 

4 


1078 

4 


4 

4 

4 

FA0 


C8 


7  DC 
190 
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CDC  IBM 
(Octal)  (Hex) 


S¥MEM(50,4)  The  output  characteristics  of  SWMEM  310  640 

nonlinearities.  S WMi=S WMEM (i,  3) 

-SWMEM(i,  2)*(SWMEM(i,  1)- “input"). 

SWMEMfi,  4)  is  a  flag  indicating  the 
state  and  a  change  of  state  in  SWMEMi. 

MAXSWS  The  maximum  i  for  which  SWCHi  is  14 

serviced 

MAXMEM  The  maximum  i  for  which  SWMi  is  14 

se  rviced 

NEVENT  The  number  of  EVENTs  to  be  serviced,  1  4 

set  by  the  user  on  the  *NEVENT  card. 

C0MM0N/UNIP1  /  VPL0T ( 100) ,  N0PL0T,  NP0INT,  NL0CAL, 

MAX  150  330 

VPL0T(1OO)  Temporary  storage  for  the  current  144  320 

(100)  PL0T  variables 

N0PL0T  A  flag  to  prevent  saving  of  any  PL0T  1  4 

variables.  If  N0PL0T*O,  no  PL0T 
variables  are  saved  and  0UTPUT  is  only 
called  at  print  times.  (N0PL0T  =  O 
nominally) 

NP0INT  During  run  time,  the  actual  number  of  14 

points  saved  on  TAPE  11  for  plotting. 

At  end  of  run,  the  act\’"l  number  of 
points  to  be  plotted. 

NL0CAL  During  run  time,  the  number  of  points  in  1  4 

the  PL0T  buffer.  (NL0CAL  <  2000  for 
CDC,  4000  for  IBM.  )  At  conclusion  of  run, 
the  number  of  frames  per  plot  point.  Must 
not  be  changed  by  user. 

MAX  The  number  of  words  per  plot  frame  1  4 

written  onto  TAPE  11  by  ESP 


C0MM0N/UNIP2/HEAD(6O),  0UT(6O),  TITLE  (8).  TODAY, 
NHEAD,  LINES,  NPA GE  ,  J  LINE,  MX  L,  IF0RM(3) 


HEAD(60) 


0UT(6O) 


TIT  LE  (8) 


T0DAY 


NHEAD 


LINES 


NPAGE 


J  LINE 


MXL 


IF0RM(3) 


Vector  which  contains  print  headings 
normally  picked  up  from  SPRINT  state- 
ment,  but  may  be  set  directly  by  over¬ 
written  F0RTRAN  statements,  or  on 
-HEADINGS  card. 

Vector  containing  output  values  to  be 
printed.  Equivalent  to  PRINT (60)  in 
0UTPUT. 

Vector  containing  title  specified  by  user 
on  -TITLE  card 

Contains  actual  date  returned  by  sub¬ 
routine  DATE  and  printed  on  output. 

The  number  of  print  variables  (head¬ 
ings) 

The  number  of  print  lines  per  block  of 
print 

Page  number  for  printout 

Used  internally  by  UNIP2  to  control 
printed  output 

Used  internally  by  UNIP2  to  control  printed 
output:  number  of  blocks  of  printout  per 
page. 

Contains  output  format  to  be  used  for 
printed  output,  based  on  accuracy  require¬ 
ments. 
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Block  Length 


CDC  IBM 
(Octal)  (Hex) 


C-4.  RESERVED  SUBROUTINE  NAMES 


NOTE 

The  subroutines  Listed  below  are  loaded  and  used  during  execu¬ 
tion  of  an  ESP  job.  The  user  should  be  careful  not  to  duplicate 
any  of  these  names  when  adding  his  own  subroutines,  except  in 
the  case  of  DERIVS,  ICC0MP,  OUTPUT,  SWINPT  SWMEMN  or 
MAIN  when  he  intends  to  supply  the  entire  routine  himself. 


ROUTINES  USED  BY  ESP  AND  GRAPH 


ABORT 

FRAMES 

L0GGRD 

SCA  LEPR 

ADAMS 

FRAMXX 

MAIN 

SECNZR 

ADMNTP 

GENGRD 

NABLE 

SHIFT 

AND 

GRAPH 

NEXTCHR 

SKIPFIL 

BUFF 

GRAPH2 

NEWGRD 

SKPFIL 

CKBLNK 

GRAPHX 

NOTIFY 

STDGRD 

COMPL 

ICC0MP 

NUMBER 

SWINIT 

CONS 

ICKBLNK 

NUPLOT 

SWINPT 

DEC0D 

IDEC0D 

NX  T  C  HR 

SWMEMN 

DERIVS 

IDFRAM 

0R 

SWTCHE 

E0FS1M 

IPICK 

Output 

SYMBOL 

ENC0D 

JUNK 

PACKER 

SYSTEM  = 

EOF 

LABCHK 

PARRAY 

TIM2G0 

ESPCTL 

JL1BRST 

PINOUT 

TIME  IN 

ESPII 

L1BSET 

PLOTS 

TIMEOU 

ESPRNT 

L1ERR 

PLTSYM 

TIMEOUT 

ESPLOT 

LEVEL 

READIT 

EVENTS 

LEVEL1 

REMARK 

FILBUF 

LEVEL2 

REST0R 

FILLBUF 

LINGRD 

SCA  LEP 

C-ll 


ROUTINES  USED  BY  WHELP 


*: 

Pi 

t 


y- 


CROSS 

MATINV 

MATZR0 

SCAMAT 

IDENT 

MATMAT 

MOVE 

TRNSML 

MATADD 

MATSUB 

NEGATE 

TRNSPS 

COMMON 

BLOCK  NAMES  (may  not  be  used  as 

subroutine  names  on  IBM) 

BASIC 

HMAXMN 

READIN 

SWHPAR 

BLANK 

LI BSC R 

RKCONT 

SWT CHS 

C0NSTS 

MISCEL 

STFPAR 

TEMSTR 

EOFSIM 

NDISPR 

STPC0N 

UNIP1 

GRAPHP 

PARS 

SWDBUG 

UNIP2 
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APPENDIX  D 


PROGRAM  CONTROL  AND  EXECUTION 

D-l.  INTRODUCTION 

How  an  ESP  program  works  can  be  considered  on  two  levels.  First, 
there  is  the  manner  in  which  the  control  cards  put  the  program  together  from 
the  user's  deck  and  the  ESP  files.  Then,  there  is  the  manner  in  which  the 
program  actually  executes  to  solve  the  user's  problem.  This  appendix  will 
attempt  to  clarify  both,  first  by  providing  a  diagram  showing  the  relationship 
of  control  cards,  compilers,  libraries  and  files  and  second  by  providing 
descriptions  and  flowcharts  of  the  major  subroutines  which  make  up  the  ESP 
library. 
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D-3.  PREC0MP 

D-3-a.  What  PRECQMP  Does 

PROGRAM  PREC0MP  is  a  precompiler,  written  in  F0RTRAN  on 
CDC  and  PLI  on  IBM,  which  reads  the  user's  ESP  language  input  deck  or  file 
and  translates  it  into  executable  F0RTRAN  routines^  to  be  used  by  the  ESP 
run-time  package.  Its  chief  functions  are  to  "crack"  the  ^control  cards 
such  as  *BL0CK,  *SWTCH,  and  *PRINT,  and  to  write  the  additional  cards 
needed  to  complete  those  subroutines  based  on  user  coding,  namely  MAIN, 
DERIVS,  ICC0MP,  OUTPUT,  SWINPT,  and  SWMEMN  (see  Appendix  D-3-b). 

PREC0MP  operates  by  making  repeated  calls  to  SUBR0UTINE 
READIT  (which  reads  the  user's  card  images)  and  by  writing  these  card 
images  out  onto  file  IMS0RC  until  a  signal  card  is  detected.  It  then  tests  the 
signal  card  to  determine  its  next  action,  which  may  be  copying  more  cards, 
setting  flags,  translating  the  data  on  the  signal  card,  or  writing  additional 
F0RTRAN  statements  onto  IMS0RC  to  complete  the  subroutines.  Since 
READIT  depends  upon  blanks,  $  terminators,  and  *  to  delimit  fields  of  data 
and  to  tell  it  how  to  handle  data,  formats  for  all  ESP  cards  should  be  followed 
carefully. 

Each  time  PREC0MP  is  executed,  it  will  do  the  following  ten  steps 
in  order,  although  substeps  may  be  in  any  order,  as  indicated: 

1.  Call  TIMEIN  to  get  precompiler  starting  time. 

2.  Process  ^METHOD  card  if  used. 

3.  Write  PR0GRAM  MAIN,  specifying  proper  integration 
package,  or_  copy  user's  PR0GRAM  MAIN,  if  provided. 

4.  Copy  all  user- supplied  subroutines  and/or  functions. 

— 

If  WHELP  statements  are  used,  however,  they  are  copied  as  is  and  must  be 
converted  to  FORTRAN  by  the  WHELP  precompiler. 
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5. 


Write  SUBROUTINES  DERIVS,  ICCOMP,  and  OUTPUT  in 
any  order  as  follows: 


a.  Write  SUBROUTINE  ICCOMP  using  all  cards  contained 
between  *ICCOMP  and  ’-ENDIC.  If  no  ’-ICC0MP  is 
used,  write  a  dummy  routine, 

b.  Write  SUBROUTINE  OUTPUT  using  all  cards  contained 
between  *OUTPUT  and  *END0UT.  If  no  ^OUTPUT 
appears,  use  data  on  SPRINT  card.  If  neither  ’-OUTPUT 
nor  SPRINT  is  used,  write  a  dummy  routine. 

c.  Write  SUBROUTINE  DERIVS  using  all  cards  contained 
between  *  DERIVS  and  ’-ENDDERIVS.  Within  this 
section,  do  the  following  in  any  order: 

(1)  Copy  FORTRAN  and  WHELP  statements  as  given. 

(2)  Translate  ’-BL0CK  cards  into  FORTRAN  and 
write  as  part  of  DERIVS. 

(3)  Process  ’-SWTCH  and  ’-SWMEM  cards  by  writing 
SUBROUTINES  SWINPT  and  SWMEMN  containing 
the  input  expressions  and  by  adding  code  to 
DERIVS  to  define  SWTCH  and  SWMEM  output. 

6.  Write  ’-HEADINGS,  ’-SWTCHES,  and  *SWMEMCNT  cards  on 
TAPE  12. 

7.  Copy  remaining  input,  such  as  ’-IV,  ’-PAR,  or  ’-RUN  cards 
and  user  data  cards,  up  to  EOF,  onto  TAPE12. 

8.  ENDFILE  12  and  REWIND. 

9.  ENDFILE  IMSORC  and  REWIND. 

10.  Call  TIMEOUT  to  compute  precompiler  time  used. 

D-3-b.  CARDS  WRITTEN  BY  PREC0MP 


Below  is  a  listing  of  the  cards  written  by  PRECOMP  which  are 
added  to  the  various  segments  of  the  user's  coding  to  produce  complete  sub¬ 
routines.  If  the  user  chooses  to  provide  any  or  all  of  these  routines  himself, 
he  should  be  careful  to  include  all  of  the  cards  listed  and  to  insert  his  own 
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coding  where  indicated.  Notice  that  when  a  particular  routine  has  no  function 
and  is  written  as  a  dummy  (see  SUBR0UTINE  ICC0MP  in  Example  Problem, 
Section  II)  not  all  of  the  cards  below  will  be  listed.  Also,  if  the  WHELP 
precompiler  is  used,  it  writes  additional  common  block  statements  for  its 
own  needs . 


MAIN  PROGRAM 


PR0GRAM  MAIN  (TAPEll,  TAPE  12,  INPUT  =  TAPE  12,  0UTPUT) 
EXTERNAL  DERIVS,  METH0D,  INTERP 
CALLESPII  (DERIVS,  METH0D,  INTERP) 

END 

where 


METH0D= 

INTERP= 


Adams 

ADAMS 

ADMNTP 


RK2 

ESPRK2 

INRKPC 


RK4  Predictor /Cor  rector 
ESPRK4  ESPPC 

INRKPC  INRKPC 


DERIVATIVE  SUBROUTINE 


SUBR0UTINE  DERIVS (T,  Y,  DY,  ST0P) 

DIME NSI0N  Y(100),  DY(100),  PAR(IOO) 

C0MM0N/SWT CHS / SWT CH (50) ,  SWMEM(50,  4),  MAXSWS,  MAXMEM,  NEVENT 
C0MM0N/ PARS /PAR 

[coding  which  defines  derivative  equations  as  DY's  and  the  desired 
outputs  of  any  discontinuities  used. 

RETURN 

END 


*SWTCH  INPUTS  SUBROUTINE 


SUBR0UTINE  SWINPT  (VALUES,  T,  Y) 

DIMENSI0N  VALUES(l),  Y(l),  PAR(IOO) 

C0MM0N/SWTCHS/SWTCH(5O),  SWMEM(50,  4),  MAXSWS,  MAXMEM,  NEVENT 
C0MM0N/PARS/PAR 

[coding  which  defines  input  expressions  to  SWTCHs  and  stores  them 
in  array  VALUES. 

RETURN 

END 


D-6 


»SWMEM  INPUTS  SUBROUTINE 

SUBROUTINE  SWMEMN  (VALUES,  T,  Y) 

DIMENSION  VALUES(l),  Y(I),  PAR(IOO) 

COMMON/SWTCHS/SWTCHCSO),  SWMEM(50,  4),  MAXSWS,  MAXMEM,  NEVENT 
cOmmOn/pars/par 

[coding  which  defines  input  expressions  to  SWMEMs  and  stores  them 
in  array  VALUES. 

RETURN 

END 


OUTPUT  SUBROUTINE 


subroutine  Output  (t,  y,  dy,  plot,  print, stop) 

DIMENSION  Y(100),  PAR(IOO),  PLOT(10),  PRINT(60),  DY(100) 
COMMON/SWTCHS/SWTCH(50),  SWMEM(50, 4),  MAXSWS,  MAXMEM,  NEVENT 
COMMON/PARS/PAR 

["coding  which  defines  print  and  plot  values  and  stores  them  in 
[_PRINT  and  PLOT,  respectively. 

RETURN 

END 


INITIAL  COMPUTATIONS  SUBROUTINE 


SUBROUTINE  ICCOMPfT,  Y) 

DIMENSION  Y  (100),  PAR(IOO) 

COMMON/SWTCHS/SWTCH(50),  SWMEM(50,  4),  MAXSWS,  MAXMEM,  NEVENT 
COMMON/PARS/PAR 

[coding  which  defines  any  initial  conditions,  computes  program 
constants,  or  performs  any  task  involved  with  program  initialization. 
RETURN 
END 


D-4.  RUN-TIME  ROUTINES 

Once  PRECOMP  is  finished,  IMSORC  will  contain  PROGRAM  MAIN, 
SUBROUTINES  ICCOMP,  OUTPUT,  DERIVS,  SWINPT,  SWMEMN,  EVENTS 
and  NOTIFY,  and  any  other  subroutines  provided  by  the  user  to  his  program. 

In  order  to  complete  the  program  and  make  it  executable,  a  group  of  sub¬ 
routines  to  be  referred  to  as  "run-time'1  routines  will  be  selected  from  the 
ESP  library  and  added  to  the  program.  The  internal  workings  of  most  of  these 
routines  probably  are  not  relevant  to  the  user,  but  a  brief  description  of  each 
follows. 
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To  aid  the  user  in  understanding  ESP  and  perhaps  in  debugging  his 
program,  some  further  information  is  included.  Appendix  D-4-b  shows  the 
overall  relationship  of  subroutines  during  execution.  Further,  Appendices 
D-4-c,  i-vi,  contain  schematic  flowcharts  of  those  run-time  routines  most 
significant  in  program  control  and  logic,  namely,  SUBROUTINES  ESPII  and 
ESPCTL,  and  the  integration  routines  ADAMS,  ESPRK4  (ESPRK2),  and 
ESPPC . 

D-4-a.  Routines  Provided  by  ESP  (Run-Time  Routines) 


ESPII 


ESPCTL 


0DESOL 


INTERP 


ESPRNT 

ESPLOT 


Controls  overall  execution  by  such  operations  as  reading  and 
interpreting  run-time  cards,  controlling  multiple  cases,  and 
calling  plot  routines  (see  flowchart,  Appendix  D-4-c-ii). 

Controls  all  of  the  tasks  needed  to  execute  one  :>;RUN  card, 
which  includes  initializing  and  printing  variables,  calling  the 
integrator  routine  selected,  printing  warnings  if  integration 
was  unsuccessful,  checking  for  switches  and  calling  the  appro¬ 
priate  switch  routines,  and  storing  print  and  plot  data  at  the 
correct  times  (see  flowchart..  Appendix  D-4-c-iii). 


The  integration  routine,  which  will  be  one  of  the  following: 
ESPPC  Predictor-corrector  method 

ESPRK2  Second  order  Runge-Kutta 

ESPRK4  Fourth  order  Runge-Kutta 

ADAMS  Adams  integration 


(See  Chapter  IV,  Integration  Package,  and  Appendices 
D-4-c,  iv,  v,  vi.  ) 


The  interpolation  routine  used  to  interpolate  data  for  printout 
and  switchings,  which  will  be  one  of  the  following: 

INRKPC  Used  for  predictor-corrector,  RK2  or  RK4 

ADMNTP  Used  for  Adams  integration 

Calls  OUTPUT  to  obtain  print  data  and  does  actual  printing 
of  output. 

Calls  0UTPUT  to  obtain  plot  data  and  stores  plot  data  for 
later  plotting. 
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SWINIT 

SWTCHE 

SECNZR 

DATE 

READIT 

TIMEIN 

TIM0UT 

GRAPH 

RESTQJR 

FILBUF  ) 
PACKER  } 
,  SKPFIL  | 

DIFTAB 


Initializes  switches  and  reinitializes  them  after  a  switching. 

Evaluates  SWTCH,  SWMEM,  and  EVENTS  inputs  by  calling 
SWINPT,  SWMEMN,  and  EVENTS,  detects  sign  or  state 
changes,  locates  zero  crossings,  and  flags  outputs. 

Finds  the  zero  crossing  when  SWTCHS,  SWMEMs,  or 
EVENTS  have  been  detected. 

Returns  date  on  which  run  is  executed. 

Reads  data  from  specified  input  stream,  terminating  on  the 
character  indicated  (blank  or  $). 

Records  solution  starting  time. 

Records  solution  stop  time. 

Does  actual  plotting. 

Used  to  manipulate  plot  buffers. 


Used  internally  for  file  reading  and  manipulation. 


Called  by  ESPPC  to  see  if  stepsize  doubling  will  introduce 
numerical  instability. 
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D-4-b. 


(  MAIN 


Relationship  of  Routines  During  Execution 


NOTES: 


K  ESPII  > 


K  skpfil) 
»(~timein) 

ESPCTL  ~)~ 


h<™7>i 


— DATE 
<  IEADIt) 


K  timout) 


►TrestorJ 

j-*C  EILBUF  ) 


ICCOMP 


K 


EVENTS 


DERIVS 


|— »(  packer) 


GRAPH  } 


1.  Q  NAME  ^  -  routine  is  part  of  ESP 
j  NAME  |  - 


SWINPT 


SWMEMN 


-c  FSPL0T> 

— ESPRNT^- 


OUTPUT 


OUTPUT 


►(oDESQl)- 


DERIVS 


— »fsWTCHE)- 
— c  NTERP  ) 


I— »C  READ  IT  ) 


K  I NTERP  ) 


swinit)— 


SWMENN 


EVENTS 


NOTIFY 


SWINPT 

p 

SWMEMN 

r*(  1  NTERP  ] 

- 

EVENTS 

- 

SWMEMN 

SWINPT 

:sprnt 

K 


OUTPUT 


ESPLOT 


OUTPUT 


routine  contains  user  code 


2,  Many  routines  are  called  more  than  once  by  their 
calling  routine  but  are  shown  only  once 

3.  ODESOL  and  I  NTERP :  are  internal  variables  representing  the  following 
subroutine  names,  depending  on  the  integration  method  selected: 


Method 

ODESOL 

1 NTERP 

Predictor /Corrector 

ESP  PC 

INRKPC 

RK2 

ESPRK2 

INRKPC 

RK4 

ESPRK4 

INRKPC 

ADAMS 

ADAMS 

ADMNTP 
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D-4-c.  Flowcharts  of  Significant  Routines 

D-4-c-i.  Explanation  of  Flowchart  Conventions 


»(  OUTPUT  } 


Call  the  subroutine  whose 
FORTRAN  name  is  OUTPUT 


Store  plot  data 


Call  subroutine  UNI  PI,  which  performs 
the  function  described  in  the  box  beneath 
it.  The  two  way  arrow  indicates  that 
program  control  returns  to  the  main  line 
at  completion  of  the  subroutine 


Compute  H  only  if  H  <  0;  then 
continue  on  "NO"  path 


1 


Enter  at  this  point  after  a  branch 
from  another  location 


Branch  from  this  point  to  the 
corresponding  entry 


l  P-4 
H  ■  0. 5  *  H  I 


Multipage  flowcharts  continue  from  the 
dangling  arrow  in  the  lower  left  corner 
of  a  page  to  the  entry  arrow  at  the 
upper  left  of  the  next  page 


FIRSTP 


Names  appearing  in  all  capital  letters 
represent  actual  FORTRAN  names 


D-ll 


D-4-c-ii 


Flowchart  of  Subroutine  ESPII 
(page  1  of  1 ) 


D-4  -c  -iii. 


Flowchart  of  Subroutine  ESPCT  L 
(page  1  of  3) 
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D-4-c-iv. 


Flowchart  of  Subroutine  ADAMS 
(page  1  of  2) 


KOLD  ~ JSTART 
CRASH  =  KFLAG 


D-l 


6 


8 


D-4-c-v.  Flowchart  of  Subroutine  ESPRK4  (ESPRK2) 
(page  2  of  2) 
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1 


D-4-c-vi,  Flowchart  of  Subroutine  ESPPC 
(page  1  of  3) 


Do  single  interval 
P/C  integration 


JSTART  ■  1 

nz 


RETURN 


►C  PER  I  VS  ) 


NODUBL  *  I 

Compute  Derivs 

I  ! 

at  T  +  H 

ESPRK4  J 


Take  2  steps  using  RK4 


.  YES 

Take  1  unchecked 

>  JL 

step  with  RK4 

NUMSTP  *  NUMSTP  +  1 
ICOUNT  *  -3 
KCOUNT  •  0 
DOUBLE  *  .FALSE. 


DOUBLE  *  .FALSE. 


Do  double  interval 
P/C  integration 


Compute  Derivs  at  T  +  2  H  | 


YES 

(Stay  with  single  interval  1 

LCOUNT  *  -4 

Double  the  interval 
L  COUNT  =  -2 
KCOUNT  -  0 
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APPENDIX  E 


INTEGRATOR  EQUATIONS 


In  all  descriptions  given  below,  assume  a  differential  equation  of  the 

form 

=  f(t’y) 

with  v  a  vector.  For  ease  of  notation,  y  =  y(t  ). 

E-l.  ADAMS  INTEGRATION 

Adams  integration  is  a  highly  complex  variable -order ,  variable  -  step 
algorithm,  which  is  completely  documented  in  Ref.  6.  As  its  complexity 
precludes  condensation,  the  user  is  referred  to  Section  IV-B  and  Appendix 
D-4-c-iv  for  an  overview  of  the  method  and  to  Ref.  6  for  the  actual  algorithm. 

E-2,  SECOND-ORDER  RUNGE-KUTTA 

yn  +  l  =  yn  +  °’5h(k0  +  kl}  (E’1} 


where 

k0  =  f(t„’  V 

kl  *  f(t„  +h’yn  +hk0> 

In  the  RK2  fixed-step  mode,  Eq.  (E-l)  is  used  to  take  two  steps  at  a 

time  before  checking  for  such  items  as  print,  plot,  and  switchings.  In  the 

RK2  variable  -  step  mode,  a  step  of  size  2h  is  taken  first  and  compared  with 

the  result  of  two  normal  steps.  If  y  is  the  result  of  the  2h  step  and  y  _ 

r  "n+2  1  n+d 

is  the  result  of  2  normal  steps,  the  estimated  error  vector  is 

err  =  <yn+2  ‘  yn+2)/3*0 


E-l 


which  is  added  to  V n+2  to  '•mProve  the  accuracy.  This  error  vector  is  used 
to  control  stepsize  based  on  the  test  outlined  in  Section  IV-C,  namely,  let 

bnd  =  eps  X  max  (y^+2>  q) 

then 

1.  If  err  >  bnd  in  any  component,  halve  the  stepsize  (if  allowed) 
and  retry. 

2.  If  err  <  bnd/30.  0  in  every  component,  set  h  =  min(stepmax,  2h) 
for  the  next  step. 

E-3.  FOURTH-ORDER  RUNGE-KUTTA 

yn  +  i  =  yn  +  h(k0  +  Zkl  +  2k2  +  k3)/6.  0  (E-2) 

where 

kn  *  f(t  ,  y  ) 

0  n’  ’  n 

k.  =  f(t  +  0.  5h,  y  +0.  5hk„) 
in  '  n  U 

k.,  s  f (t  +  0.  5h,  y  +0.  5hk. ) 

2  n  ’n  1 

kQ  =  f(t  +  h,  y  +  hkQ ) 

3  n  Jn  2' 

In  the  RK4  fixed-step  mode,  Eq.  (E-2)  is  used  to  take  two  steps  at  a  time 
before  checking  for  such  items  as  print,  plot,  and  switchings.  In  the  RK4 
variable -step  mode,  a  step  of  size  2h  is  taken  first  and  compared  with  the 
result  of  two  normal  steps.  If  yn+2  is  the  result  of  the  2h  step  and  yn+2  *-s 
the  result  of  the  two  normal  steps,  the  estimated  error  vector  is 

e"  =  ^n+Z  -?„+2)/15-0 

which  is  added  to  yn+2  imProve  the  accuracy.  This  error  vector  is  used 
to  control  stepsize  based  on  the  test  outlined  in  Section  IV-C,  namely,  let 

bnd  =  eps  X  max(yn+2>  q) 


E-2 


then 


1. 


If  err  >  bnd  in  any  component,  halve  the  stepsize  (if  allowed) 
and  retry. 


2.  If  err  <  bnd/  150.  0  in  e ve ry  component  set  h  =  min  (stepmax,  2h) 
for  the  next  step. 


E - 4.  HAMMING  PREDICTOR  CORRECTOR 


Once  four  back  values  have  been  created  using  Eq.  (E-2),  the  following 
formulae  are  used  at  each  step  (primes  indicate  derivatives). 


p  ,  =  y  ,  +  4h ( 2 y 1  -  y1  .  +  2y‘  )/3 

rn  +  l  7n-3  7n  ■'n-l  7n-2 


m„xi  =  H 2 (p  -  c  )/121 


n  n 


n  +  1  *n  +  l 

Cn+1  =  [9yn  '  yn-2  +  3h(m;  +  l  +  2yn  "  yn-l)]/8 


1 


The  difference  from  the  previous  step  is  divided  by  32  to  account  for 

halving,  multiplied  by  32  to  account  for  doubling,  and  set  to  zero  following 
an  RK4  restart. 

Stepsize  doubling  is  only  attempted  if  err  <  bnd/100  and  the  number  of 
successful  predictor-corrector  steps  has  been  at  least 

1.  3  after  a  RK4  restart  or  halving 

2.  2  after  a  successful  doubling 

3.  4  after  a  doubling  failure 
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CASES:  MULTIPLE  RUNS 
AND  LARGE  SIMULATIONS 


APPENDIX  F 


SPECIAL  CASES:  MULTIPLE  RUNS  AND  LARGE  SIMULATIONS 


F-l.  MULTIPLE  RUNS 

F-l-a,  Multiple  Runs  Varying  Run-Time  Data  Cards 

If  a  series  of  job  runs  is  to  be  made  in  which  the  only  changes  from 
one  run  to  the  next  are  in  items  which  can  be  input  on  run-time  data  cards, 
then  any  number  of  runs  can  be  made  as  one  job.  The  necessary  run-time 
cards  are  simply  stacked  in  sequence,  according  to  the  following  rules: 

•  IVs,  PARs,  SWMEMDATA,  Qs,  EPS,  and  all  other  run 
time  data  cards  except  SWMEMSET  retain  their  values 
until  they  are  reset  by  the  user's  program  or  on  a  new 
run-time  card  such  as  ’-IV  or  -PAR. 

•  ’-SWMEMSET,  if  used,  must  be  redefined  for  each  run 
since  it  is  changed  during  execution. 

•  Each  ’-RUN  card  produces  reexecution  of  the  program  as 
soon  as  it  is  encountered,  so  must  always  be  the  last  run¬ 
time  card  of  a  case. 

•  A  set  of  *GRAPH  cards  must  follow  each  ’-RUN  from  which 
plots  are  expected. 


EXAMPLE: 

’-ENDIC 

*IV  0.5  Y3=0. 01  $ 

-PAR  5.0  3.  14  57.6  10.0  $ 

-SWMEMDATA 

'  First  run,  with  plots. 

1  1.0  2.  0  0.7  0.  7  1.  0  $ 

-SWMEMSET  1  $ 

-EPS  A  LL=1 . 0E-  10  $ 

-RUN  3  0.  1.  0  100.  $ 

’-GRAPH  1  3 

Second  run:  reset  SWMEM  in 

*PAR  10.0  P4=20.  0  $ 

saturation;  change  PARs,  all 

-'SWMEMSET  1  $ 

others  the  same;  no  plots. 

’-RUN  3  0.0  1.0  100.  $ 

’-IV  0.1  0.1  0.1  $ 

’-SWMEMSET  1  $ 

*RUN  3  0.  1.0  100.  $ 

*GRAPH  2  3 

• 

Third  run:  reset  SWMEM  in 
saturation,  change  IVs,  use 
PARs  from  second  run,  all 
others  from  first,  make  plots. 
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F-l-b. 


Multiple  Runs  with  the  Same  Run-Time  Data  Cards 


Sometimes  the  only  changes  from  run  to  run  are  in  the  data  read 
from  the  user's  file  (see  case  below)  and  the  user  has  a  rather  lengthy  list  of 
run-time  cards  (for  example,  lengthy  *GRAPH  cards)  which  he  prefers  not  to 
duplicate  for  each  of  the  stacked  cases.  This  can  be  avoided  in  the  following 


manner: 


F-l-c. 


•  Use  a  ^RETURN  card  as  the  last  run-time  card,  after  *RUN 
and  all  *GRAPH  cards. 

•  Have  ICC0MP  read  the  data  from  the  user's  file  and  terminate 
program  execution  when  all  data  is  exhausted. 

•  Write  the  MAIN  program,  being  sure  to  declare  the  user's 
file  on  the  PROGRAM  card,  and  include  logic  to  backspace 
TAPE12  (this  is  the  file  on  which  run-time  cards  reside 
during  execution)  after  each  case  exactly  as  many  command 
cards  as  are  needed  for  one  case  and  then  loop  back  to  the 
call  to  ESP.  (If  all  *  command  cards  are  used  for  each  case, 
TAPE12  may  be  rewound  instead  of  backspaced. ) 

Multiple  Runs  Varying  Data  to  be  Read  In  by  User 

If  the  user  has  data  decks  he  wishes  to  read  in  from  ICC0MP  and  he 


wants  to  stack  up  a  number  of  cases,  he  may  use  the  same  stacking  of  run¬ 
time  cards  as  shown  in  the  above  example,  but  must  also  do  several  other 
things : 

•  Either  stack  his  data  cards  for  each  case  immediately 
following  the  *RUN  to  which  they  apply  (see  Section  VII-E) 
or  place  all  of  the  data  on  a  separate  file  before  running  the 
ESP  job  and  declare  this  file  on  the  PR0GRAM  MAIN  card 
by  writing  his  own  PR0GRAM  MAIN. 

•  Read  the  data  in  from  ICC0MP  by  means  of  READ  or 
NAMELIST  statements.  [If  he  plans  to  test  for  the  end  of 
data  for  a  given  case  by  using  IF  (E0F.  .  . ,  he  must  be  sure 
to  write  E0Fs  on  his  data  file  when  he  creates  it,  by  using 
the  F0RTRAN  ENDFILE  n  between  data  sets.] 


F-l-d. 


Making  Plot  Overlays  of  Data  from  Multiple  Runs 


Normally  the  tape  containing  the  plot  data  is  always  rewound 
whenever  a  *RUN  card  occurs.  However,  by  using  *RUNC  in  place  of  *RUN 
for  cases  after  the  first,  the  user  may  cause  subsequent  plot  information  to  be 
written  onto  the  file  following  the  plot  data  from  the  previous  case(s).  To 
retrieve  this  new  information  for  plotting  the  user  must  specify  the  appropriate 
case  number  in  parentheses  following  the  *GRAPH. 


EXAMPLE: 


-MAX  PLOTS 

25 

!RUN  12  0 

1  10 

$ 

:PAR  P86  = 

:  34.  5  $ 

:RUNC  12 

0  1  10 

$ 

=  PAR  P72  = 

:  12.5  $ 

-RUNC  12 

0  1  10 

$ 

= GRAPH  1 

2 

PRINTER  PLOT  OF  DATA  FROM  FIRST  CASE 
-GRAPH  (2)  1  2 

PRINTER  PLOT  OF  DATA  FROM  SECOND  CASE 
-GRAPH  1  3  TYPEF 

FILMPLOT  OF  DATA  FROM  CASES  1-3 
-GRAPH  (2)  1  3  TYPEF  OVERLAY 

-GRAPH  (3)  1  3  TYPEF  OVERLAY 


F-l-e.  Running  the  Solution  Backward  and  Forward 

Boundary  value  problems  and  other  problems  where  it  is  desired  to 
have  the  capability  of  running  the  independent  variable  in  either  direction 
may  be  handled  by  having  one  PAR,  say  PAR(99),  be  +1,0  for  forward  solu¬ 
tion  and  be  -1.0  for  backward  solution.  Thus  use 

T  =  PAR  (99)  *  T 

I  Expressions  defining  thel 
I  derivatives  in  DY  I 

T  =  PAR (99)  *  T 
DO  n  i=l,  neq 

n  DY (i)  =  PAR (99)  *  D Y(i) 


in  the  derivative  routine  and  an  arrangement  such  as  the  following  for  the 
run  time  cards: 


-PAR  P99  =1.0  $ 

*RUN  3  2.0  0.  1  10.0  $ 

*PAR  P99  =  -1.0  $ 

*RUN  3  -10.0  0.1  -2.0  $ 


forward  solution 
backward  solution 


It  will  also  be  necessary  to  copy  DY  into  Y0  in  ICCOMP  by  including: 

COMMON/BASIC/T0,  TF,  TP,  Y0 (100),  YPRNT(IOO),  DY(100) 
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USING  ESP  FOR  LARGE  SIMULATIONS 


F-2-a.  Maximum  Dimensions 


In  general,  any  combination  of  ESP  variables  and  special  facilities 
may  be  used  within  one  program.  However,  each  of  the  following  items  is 
limited  to  the  total  number  indicated: 


Derivatives: 


Discontinuities: 


Print  Variables: 


Plot  Variables: 


Parameters : 


100,  whether  defined  as  DYs  or  in 
*BL0CK  form. 

1 50  total 

50  defined  as  *SWTCH 
50  defined  as  *SWMEM 
50  defined  as  EVENTS 

60  defined  by  *PRINT  plus 

any  number  of  variables  that  are  user- 

formatted. 

100  if  using  PL0T  and  :':GRAPH  or 
any  number  if  user  writes  his  own  plot 
file  and  uses  other  means  of  plotting. 

100  stored  and  passed  by  PAR  array 
plus  any  number  sto’ed  and  passed  by 
user. 


F-2-b.  Maintaining  Flexibility 

Since  especially  large  simulations  often  require  changes  and  revi¬ 
sions,  it  is  highly  desirable  to  structure  them  in  a  manner  which  makes 
additions  and  deletions  as  painless  as  possible.  Below  is  a  suggestion  for 
one  method  of  maintaining  flexibility  in  numbering  and  referencing  the  deriva¬ 
tives,  which  the  user  may  find  adaptable  to  his  program. 

The  basic  goal  is  to  start  with  a  structure  which  eases  the  problem 
associated  with  the  inevitable  changes  required.  The  approach  suggested 
here  is  to  modularize  and  to  use  pointers  such  that  the  modules  have  maxi¬ 
mum  independence. 


•  Define  a  common  block  containing  two  arrays,  each  having 
at  least  as  many  words  as  there  are  modules,  e.g., 
C0MM0N/IP0INT/IP0INT  (50),  NL0CAL(5O) 

•  Set  NL0CAL(I)  equal  to  the  number  of  derivatives  defined 
in  the  Ith  block. 

•  Set  IP0INT(1)  =  0  and  define  the  remainder  of  the  IP0INTs  by 

IP0INT  (I)  =  IP0INT  (1-1)  +  NL0CAL  (I-  1 ) 
for  1=2, .  .  . 

•  Within  the  Ith  module  (it  need  not  be  a  separate  subroutine) 
define  the  DYs  by,  say, 

L0C  =  IP0INT  (I) 

DY  (L0C  +1)  =  ... 

DY  (L0C  +2)  =  ... 

• 

With  this  scheme  one  need  only  know  the  correspondence  of  physical 
variables  within  a  module.  Thus,  if  the  angular  displacements  of  body  5  were 
the  fourth,  fifth,  and  sixth  variables  within  module  5,  they  could  be  used 
anywhere  else  by 

L0C  =  IP0INT  (5) 

ADI  =  Y  (L0C  +  4) 

AD2  =  Y  (L0 C  +  5) 

AD3  =  Y  (L0C  +  6) 

This  approach  may  also  be  useful  if  the  number  of  Y's  varies  with 
the  input  such  as  is  encountered  in  structures  programs. 

F - 2 - c .  Producti on  Runs  with  a  Compiled  Program 

If  produc.ion  runs  are  to  be  made  with  an  ESP  program  which 
requires  considerable  time  for  PREC0MP,  WHELP,  and/or  FTN  compilation, 
it  may  be  desirable  to  create  a  binary  version  of  the  program  and  a  separate 
TAPE  12  file  so  that  runs  can  be  made  without  recompilation.  This  may  be 
done  either  with  cards  or  files  via  the  following  steps: 


F-b 


•  Using  the  usual  control  cards  for  an  ESP  (and  WHELP) 
program,  compile  the  source  program  and  either  punch 
out  the  LGO  file  (binary)  or  catalog  it  as  a  permanent  file. 

•  Add  to  the  run-time  data  cards  those  cards  normally  written 
by  PREC0MP,  namely,  *SWTCHES,  *S WMEMCNT,  and 
-^HEADINGS,  so  that  the  run-time  data  section  looks  like  this 
(These  cards  also  may  be  used  as  a  deck  or  put  on  a  perma¬ 
nent  file.  ): 

❖SWTCHES  n  [where  n  is  actual  number 

of  *SWTCHs  used.] 

:'-SWMEMCNT  n  [where  n  is  actual  number 

of  *SWMEMs  used.] 

^HEADINGS  name  ....  name  $  [where  name ....  name 

1  m  .  .  1  .  ,  m 

are  actual  print  headings 

specified  in  *PRINT  state¬ 
ment.  If  SPRINT  is  not 
used,  this  card  is  unneces¬ 
sary.  ] 

(Any  optional  run-time  cards,  such  as  *PAR,  *IV,  etc.  ) 

*RUN.  .  .  [usual  format] 

*GRAPH..  .  [usual  format] 

♦RETURN 
To  make  runs: 

•  If  the  compiled  program  and  run-time  cards  are  on  files, 
simply  attach  the  program  file  and  call  it  L G<0,  attach  the 
TAPE12  file  and  call  it  TAPE1Z,  and  then  execute  LG0). 
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Example : 

ATTACH(LIB1,  2NEWRESP) 

ATTACH(LIB2,  3FTNPL0TUB) 

LIBRARY  (LIBl,  LIB2) 

ATTACH(LG0,  2BINAR  YPRG,  ID  =  VAP185) 
ATTACH(TAPE12,  2 TAPE  12,  ID=VAP185) 

_ LG0. _ 

•  Alternatively,  the  binary  card  deck  and  the  run-time  cards 
may  be  used  in  the  following  setup: 

ATTACH(LIB1,  2NEWRESP) 

ATTACH  (LIB 2, 3FTNPL0TLIB) 

LIBRARY  (LIB1,  LIB2) 

COPYS (INPUT ,  LG0) 

CCiPYS (INPUT,  TAPE12) 

REWIND(TAPE12) 

LG0. 

7-8-9 

[Binary  deck  of  user's  program] 

7-8-9 

[Run-time  cards  as  shown  above] 

6-7-8-9 
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DEBUGGING  SUGGESTIONS 


G-l.  TERMINATION  DURING  PREC0MP 


Check  control  cards. 

Check  deck  structure. 

Look  for  diagnostic  message  at  end  of  listing. 

Check  card  formats  for  such  items  as  required  blanks  and  $s. 

If  Precompiler  reads  off  END-QF- FILE,  check  for  $  terminators, 
especially  on  *PRINT  card. 

G-2.  TERMINATION  NEAR  BEGINNING  OF  EXECUTION 

Does  NEQ  on  run  card  agree  with  highest  numbered  derivative? 

Is  every  DY(i)  (i  <  NEQ)  defined,  even  if  just  set  =  0? 

Check  derivative  equations  carefully. 

Check  stepsize  and  any  stepsize  limits  you  may  be  specifying. 

[it  may  be  helpful  to  print  out  the  stepsize  (H).  ] 

If  system  seems  to  immediately  go  unstable,  check  HMIN  (see 
Appendix  G-6). 

G-3.  TIME  LIMIT  DURING  STARTING  PROCEDURE 


Try  running  the  program  using  all  Runge-Kutta  integration. 

Reduce  the  print  interval  to  obtain  more  printout. 

Print  the  stepsize,  H,  in  order  to  monitor  its  behavior. 

Set  HMIN  >  0.,  which  causes  integration  to  continue  by  accepting 
Y(i)'s  in  spite  of  errors. 

Check  equations  for  a  very  small  or  0.  value  of  DY  coupled 
with  a  Y0  which  is  also  very  small  or  0.  In  this  situation  the 
error  constraints  may  be  unduly  difficult  to  satisfy,  and  increasing 
the  corresponding  Q (i )  may  alleviate  the  problem. 

G-4.  EXECUTION  OCCURS  BUT  PRINTOUT  IS  ZERO  OR 
INACCURATE 


Check  PARs,  IVs,  and  all  other  inputs  as  they  are  printed 
initially. 


Are  you  attempting  to  pass  time -dependent  variables  as  PARs? 

Have  you  defined  all  output  variables  in  OUTPUT,  either  by 
computing  them  there  or  passing  them  through  common  blocks? 

Check  each  subroutine  to  make  sure  referenced  variables  are 
defined. 

G-5.  DISCONTINUITIES  ARE  NOT  WORKING  PROPERLY 


Have  you  introduced  discontinuities  only  as  SWTCHs  or  SWMEMs 
or  by  using  FIRST P  flag? 

Are  switch  inputs  properly  defined?  See  section  on  discontinuities 
to  review  restrictions. 

Are  the  allowable  timing  errors  (HSW,  HSWM,  HSWE)  suitable 
to  your  problem? 

Do  you  have  switches  driving  switches  directly?  (Review  Section  V-E) 
G-6.  "ILL-CONDITIONED  SYSTEM"  MESSAGE* 


This  generally  results  from  one  of  the  following: 

Errors  in  derivative  equations 

Discontinuities  occurring  that  are  not  at  switch  times 

Smooth  functions  which  suddenly  change  quickly  with  no  switches 
occurring 

Suggestions  for  resolving  this  problem  are  as  follows; 

Check  derivative  equations  for  coding  errors. 

Check  allowable  errors  (EPS  and  Q)  to  see  if  they  are  realistic 
for  your  problem. 

Consider  adjusting  HMIN  (default  =  0):  An  HMIN  >  0  will  cause 
acceptance  of  Y's  in  spite  of  errors  whenever  H  <  2.0*  HMIN. 
This  has  the  effect  of  forcing  integration  past  rough  spots-- 
useful  in  some  cases  but  causing  instability  in  others. 


See  Section  IV-D-1. 
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H-l.  Introduction:  How  WHELP  Works 

H-2.  Fixed  Dimension  WHELP . 

H-3.  Variable  Dimension  WHELP  .  .  . 
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WHELP 

H-l.  INTRODUCTION:  HOW  WHELP  WORKS 

WHELP  is  a  higher-level  language  preprocessor,  which  translates 
normal  (mathematical)  engineering  equations  involving  scalars,  vectors,  and 
matrices  into  F0RTRAN  statements.  It  consists  of  two  parts:  a  preprocessor 
which  converts  the  WHELP  equations  to  FORTRAN  statements,  and  a  library 
of  highly  efficient  subroutines  which  perform  the  actual  matrix  operations. 

WHELP  offers  two  big  advantages:  speed  and  accuracy.  It  saves 
coding  time  by  permitting  matrix  equations  to  be  coded  much  in  the  same 
form  as  they  are  written  (instead  of  as  a  series  of  subroutine  calls);  it  helps 
to  reduce  programmer  errors  since  the  user's  code  is  simpler  than  other¬ 
wise  required  and  maintains  a  close  resemblance  to  the  equations  it  repre¬ 
sents;  and  it  facilitates  debugging  because  any  equation  or  coding  errors  that 
do  crop  up  are  easier  to  find. 

For  example,  the  equation 


where  A,  B,  C,  D,  and  E  are 
3x3  matrices 

which  normally  the  user  would  have  to  code  as  a  series  of  subroutine  calls, 
can  be  coded  in  WHELP  with  the  one  card: 


A  =  B,  *C  +  D*E  $ 


WHELP  functions  by  first  setting  up  a  list  of  those  variable 
names  (scalars,  vectors,  and  matrices)  which  the  user  has  declared  as 
WHELP  variables.  It  then  searches  the  user's  coding  until  it  recognizes 
one  of  those  variables  on  the  left-hand  side  of  an  equal  sign.  This  is  inter¬ 
preted  as  the  signal  that  the  right-hand  side  is  a  WHELP  (vector-matrix) 
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expression  and  mast  ue  translated  to  FORTRAN.  WHELP  then  scans  the 
expression  up  to  the  $,  comparing  each  variable  to  the  list  of  declared 
WHELP  variables,  to  ascertain  which  represent  scalars,  vectors  and 
matrices,  and  to  determine  their  dimensions.  Using  the  proper  dimensions 
it  then  interprets  the  operator  symbols  (+,  *,  etc.)  to  generate  calls  to  the 
appropriate  subroutines  to  perform  matrix  multiplication,  addition,  trans¬ 
position,  etc.  The  end  result  of  the  WHELP  processor  is  an  executable 
FORTRAN  program,  in  which  the  user's  matrix  equations  will  appear  as 
comment  cards  followed  by  the  subroutine  calls  needed  to  implement  them. 
When  this  program  is  then  compiled  and  executed,  the  needed  matrix  sub¬ 
routines  will  be  loaded  from  the  WHELP  library.  The  user,  therefore,  has 
two  tasks  to  perform:  declaring  and  dimensioning  those  variables  which  will 
be  WHELP  variables  and  writing  his  vector-matrix  equations  in  a  form  that 
WHELP  can  interpret. 

Basic  WHELP  assumes  vectors  and  matrices  of  fixed  size,  but 
WHELP  may  also  be  used  with  arrays  of  variable  dimensions.  The  discus¬ 
sion  below  will  start  with  the  simplest  application,  fixed-dimension  WHELP, 
and  then  proceed  to  a  description  of  variable  dimension  WHELP  and  to  some 
special  time-saving  features  which  have  been  added  to  WHELP. 

H-2.  FIXED  DIMENSION  WHELP 

A  WHELP  variable  is  defined  to  be  any  vector  or  matrix  which 
will  be  referenced  as  a  vector  or  matrix  in  a  WHELP  expression  or  any 
scalar  which  is  to  receive  the  result  of  a  WHELP  computation.  Any  scalar 
or  single  element  of  an  array  which  will  b'  referenced  only  on  the  right-hand 
side  of  a  WHELP  expression  need  not  be  declared  as  a  WHELP  variable. 

There  are  two  types  of  WHELP  statements:  declaration  state¬ 
ments  and  WHELP  equations.  WHELP  declaration  statements  always  begin  in 
column  1  with  an  asterisk  (*)  followed  by  the  appropriate  name  (SAMESIZE, 
IDECLARE,  or  INF0RM:  See  below)  and  entries  are  terminated  by  a  dollar 
sign  ($)  preceded  by  at  least  one  blank.  WHELP  equations  always  begin  with  a 
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declared  WHELP  variable  starting  in  column  7  or  later  followed  by  an  equal 
sign  and  the  appropriate  expression  and  are  terminated  by  a  dollar  sign  ($) 
preceded  by  at  least  one  blank.  Both  may  extend  through  column  72  and  be 
continued  on  the  next  card.  No  continuation  marks  of  any  kind  are  used. 

The  statements  are  assumed  to  continue  until  the  $  terminator. 

H-2-a.  Declaring  WHELP  Variables 

Every  scalar,  matrix,  or  vector  which  is  to  be  used  as  a  WHELP 
variable  must  appear  on  a  declaration  card  within  each  (sub)  program  in  which 
it  will  be  so  used.  A  declaration  card  is  a  card  started  in  column  1  by 
’■-IDECLARE,  *SAMESIZE,  or  ^INFORM  and  terminated  by  a  $  preceded  by  at 
least  one  blank.  Alternative  forms,  which  are  both  translated  into  REAL 
statements,  are 


*IDECLARE  item  item...  item  $ 


where 


item  = 


n  (and  m) 


and 


name 
name  (n) 
name  (n,  m) 

are  integer  constants  specifying  the  number  of  rows  (and 
columns)  in  the  array 


*SAMESIZE  n  m  list  $ 


where 

n  (and  m)  are  integer  constants  specifying  the  number  of  rows  (and 

columns)  in  all  arrays  named  in  the  list.  The  parameter  n 
must  be  present,  but  m  may  be  omitted,  in  which  case  all 
arrays  named  are  singly  subscripted. 

list  is  a  list  of  names,  separated  by  at  least  one  blank,  which 

are  to  have  the  dimension(s)  given. 
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EXAMPLES: 

*SAMESIZE  3  X  Y  Z  W  $ 

*IDEC  LA  RE  F(6,  3)  G(12,  12)  D  H(6)  $ 

*SAMESIZE  12  12  A  B  C  E  $ 

These  will  be  translated  by  WHELP  into  the  FORTRAN  statements: 

REAL  X(3),  Y(3),  Z(3),  W(3) 

REAL  F(6,3),  G(12,  12),  D,  H(6) 

REAL  A(12,  12),  B(12,  12),  C(12,  12),  E(12,  12) 

Note  that  all  variables  will  be  typed  real  and  that  these  are  the  only  state¬ 
ments  needed  to  dimension  these  variables,  and  in  fact  the  variables  must 
not  be  dimensioned  elsewhere.  Also  note  that  blanks  are  the  only  delimiters 
between  list  items:  do  not  insert  commas.  Extra  blanks  may  be  inserted 
between  items. 


Note : 

1.  A  name  is  a  string  of  1-7  characters  (1-6  characters  for  IBM) 
acceptable  to  F0RTRAN  as  a  variable  name. 

2.  Embedded  blanks  are  not  allowed  in  names,  since  blanks  are 
used  as  delimiters  between  items,  but  extra  blanks  may  be 
inserted  between  items  to  improve  readability. 

3.  TEMS  and  CONSTS  are  reserved  names. 

4.  Since  both  *SAMESIZE  and  *IDECLARE  cards  are  translated 
into  REAL  statements,  they  must  precede  any  executable 
statements . 

5.  WHELP  variables  must  be  declared  in  each  routine  in  which  they 
are  to  be  used  as  such. 

6.  Scalar  variables  may  be  used  anywhere  within  WHELP  expres¬ 
sions,  but  if  they  appear  on  the  left-hand  side  of  a  WHELP 
expression  (for  example,  as  the  result  of  a  dot  product),  they 
must  be  declared  as  WHELP  variables. 
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A  third  form  of  declaration  statement,  *INF0RM,  permits 
a  submatrix  to  be  treated  as  a  WHELP  variable.  That  is,  the  variable  which 
will  appear  in  a  WHELP  expression  may  be  a  subset  of  a  vector  or  matrix, 
instead  of  the  entire  vector  or  matrix.  *INF0RM  works  much  like 
*SAMESIZE  except  no  F0RTRAN  declaration  is  written  out.  The  character 
strings  listed  on  the  *INF0RM  card  are  simply  added  to  WHELP's  list  of 
"recognized"  WHELP  variable  names,  and  the  assumption  is  made  that  these 
variables  are  dimensioned  elsewhere  by  *SAMES1ZE,  *IDECLARE,  C0MM0N, 
REAL  or  DIMENSI0N  statements.  The  format  is 


*INF0RM  n  m  list  $ 


where 


n  (and  m)  are  integer  constants  specifying  the  number  of  rows  (and 
columns)  in  all  arrays  named  in  the  list.  The  parameter 
n  must  be  present,  but  m  may  be  omitted,  in  which  case 
all  arrays  named  are  treated  as  vectors. 


list  =  iterrij  item^  item^.  .  .  item 
where 


item 


name 
name  (J) 
name  (1 ,  L) 
name  (1 ,  1,  K) 


/  The  Is  represent  a  fixed  location  for  the 
starting  point  of  the  subset.  They  could 
also  be  any  other  constant  within  the  maxi¬ 
mum  array  size  constraints,  as  long  as  the 
subset  referenced  consists  of  elements  that 
are  stored  contiguously  in  the  full  array. 


name  must  be  dimensioned  elsewhere  within  the  routine  and  the 

subscripts  J,  L,  and  K  denote  which  dimension  of  name 
is  to  be  varied  in  referencing  subsets  of  name. 


Note : 


1.  Item  is  restricted  to  10  characters  total. 

2.  Items  must  be  written  without  blanks  since  blanks  separate  items. 

3.  Only  a  totally  contiguous  subset  of  an  array  may  be  declared  a 
WHELP  variable  in  this  manner.  (See  H-3-a,  Data  Storage  and 
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Transmission.)  For  example,  the  columns  o£  a  3  x  3  matrix 
BMAT  could  be  declared  on  an  *INF0RM  card  as: 

-INFORM  3  1  BMAT(l.J)  $ 

but  not  the  rows: 

-INFORM  1  3  BMAT(J,  1)  $ 

because  the  data  in  a  matrix  row  is  not  stored  contiguously.  In 
other  words,  BMAT(J,  1)  is  the  starting  location  for  an  array  of 
the  3  next  elements  in  storage,  and  since  FORTRAN  always 
stores  a  matrix  such  as  BMAT  by  columns,  a  reference  to 
BMAT(J,  1)  where  J  =2  would  give  the  following: 


■  1 

4 

7  " 

‘2" 

if  BMAT  =  2 

5 

8 

then  BMAT  (J,  1)  = 

3 

3 

6 

9  . 

4 

EXAMPLE  (in  WHELP  code): 


C0MM0N/BLOCK1/B(5,  5) 

REAL  MM(5,  5,  3)  A(5,  5) 

*IDECLARE  M(5,  5)  X(5)  $ 

-INFORM  1  X(J)  $ 

^INFORM  5  M(1,L)  $ 

-INFORM  5  5  MM(1,  1,K)  $ 

C  STORE  THE  DOT  PRODUCT  OF  THE  JTH  COLUMN  IN  X(J). 
DO  1  J  =  l,  5 

L=J 

X(J)=M(1,  L)  .  M(l,  L)  $ 

1  CONTINUE 

C  ADD  VECTOR  X  TO  THE  LTH  COLUMN  OF  M. 

L=2 

M(l,  L)=X+M(1,  L)  $ 

C  STORE  M/K  AS  THE  KTH  SUBMATRIX  OF  MM. 

DO  2  K=l,  3 

MM(1,  1,  K)=M/K  $ 

2  CONTINUE 

USE  INFORM  TO  DECLARE  PREVIOUSLY  DIMENSIONED 
VARIABLES. 

)RM  5  5  A  B  $ 

B=M+IDENT  (A )  $ 


Note  the  following  points  in  the  example  above: 

1.  *INF0RM  1  X(J)  $  tells  WHELP  to  treat  the  string  X(J)  as  a 
WHELP  scalar.  The  1  must  appear  in  the  -"INFORM  statement. 

2.  *INF0RM  5  M(l,  L)  $  tells  WHELP  to  treat  the  string  M(l,  L) 
as  a  5-vector. 

3.  *INF0RM  5  5  MM(1,  1,  K)  $  tells  WHELP  to  treat  the  string 
MM(1,  1,  K)  as  a  5  x  5  matrix.  The  program  is  storing  1 
5x5  matrix  in  each  of  the  3  planes  of  MM. 

4.  The  maximum  dimensions  of  all  :;TNF0RM-dec  lared  variables  are 
given  elsewhere  by  DIMENSION,  REAL,  C0MM0N,  *IDECLARE, 
or  *SAMESIZE  statements. 

5.  The  *INF0RM  card  places  the  total  character  string  (for  example 
MM(1,  1,K))  in  the  table  of  recognized  variables.  Thereafter, 
the  total  string  should  be  thought  of  as  a  FORTRAN  "name"; 

that  is,  its  spelling  is  sacrosanct  and  no  changing  or  substituting 
of  variables  (for  example  K)  is  allowed,  and  it  may  be  referenced 
only  by  the  full  character  string  exactly  as  it  appears  on  the 
*  INFORM  card. 

6.  Just  as  any  WHELP  scalar  which  will  appear  on  the  left-hand 
side  of  a  WHELP  expression  must  be  declared  on  an  *IDECLARE 
card,  any  element  of  a  matrix  which  is  to  be  used  as  a  scalar  on 
the  left-hand  side  of  a  WHELP  expression  must  be  declared  on  an 
^INFORM  card  with  dimension  1,  as  in  the  example, 

♦INFORM  1  X(J)  $. 

7.  Since  "'INFORM  does  not  result  in  the  writing  of  any  dimension 
statements  (and  in  fact  will  appear  only  as  a  comment  card  in 
the  FORTRAN  listing  of  the  program),  it  may  be  used  anywhere 
within  a  (sub)program  as  long  as  it  precedes  WHELP  statement 
references  to  the  variables  it  declares. 

8.  "TNF0RM  may  be  very  conveniently  used  when  it  is  desired  to 
declare  variables  which  have  already  been  dimensioned  else¬ 
where.  (Use  of  *SAMESIZE  or  *IDECLARE  in  the  same  situation, 
since  they  result  in  REAL  statements,  would  produce  double 
dimensioning,  not  normally  allowed  in  FORTRAN.) 
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Writing  WHELP  Expressions 

A  WHELP  vector -matrix  expression  always  begins  with  a 
declared  WHELP  variable  and  ends  with  a  $,  but  otherwise  may  be  written  in 
much  the  same  form  that  it  is  written  mathematically,  simply  by  using  the 
symbols  given  in  the  table  below  for  the  desired  matrix  operations.  For 
example,  one  might  wish  to  evaluate  a  variable  after  several  coordinate  trans¬ 
formations  by 


|  XNEW }  =  a[A]  [B]  [C]  jXOLD} 

or  estimate  a  correction  term  by 

jdelxj  =  [ATA]  *  A^jr| 

Using  WHELP,  these  are  coded 

XNEW  =  ALPHA*  A*B*C*X0LD  $ 

DELX  =  INVERSE  (A,  *A)*A,  *R  $ 

and  will  appear  exactly  like  this  on  the  first  listing  of  the  program.  On  the 
FORTRAN  file  listing  of  the  program  produced  by  WHELP,  all  WHELP 
expressions  will  be  rewritten  as  comment  cards,  immediately  followed  by 
the  generated  FORTRAN  calls  (see  Calling  Sequences,  Appendix  H-2-b-ii). 

Note  that  in  the  example  above,  the  vectors  XNEW,  X0LD,  R 
and  DELX  and  the  matrices  A.  B,  and  C  must  have  been  previously  declared 
on  a  *SAMESIZE,  *IDECLARE,  or  ^INFORM  card. 
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Note  also  that  any  character  string  appearing  in  a  WHELP  expres¬ 
sion  which  is  not  identical  to  a  WHELP  declared  character  string  (ALPHA  in 
the  example  above)  will  be  treated  as  a  scalar.  That  is,  if  X0LD  is  a  declared 
WHELP  vector,  but  the  character  string  X0LD  (I)  appears  in  a  WHELP 
expression,  X0LD  (I)  is  treated  as  a  scalar.  In  general  the  acceptable  string 
length  is  1-10  characters  if  the  first  character  is  a  letter,  but  is  not  limited 
for  numbers.  However,  since  a  longer  character  string  is  sometimes 
inevitable,  for  example  ZETA  (I  +  1,  2  *  J),  WHELP  will  use  the  correct 
value  for  ZETA  (I  +  1,  2  *  J)  but  will  substitute  a  shorter  character  string 
for  the  too  long  one  when  it  writes  out  the  F0RTRAN  file.  In  any  case,  the 
total  character  string  must  not  exceed  40  characters,  or  information  will  be 
lost. 

A  WHELP  equation  is  evaluated  according  to  the  hierarchy  of  the 
operators,  given  in  Table  H-2-b-i.  In  expressions  with  like  operators, 
evaluation  occurs  from  left  to  right.  However,  as  in  standard  F0RTRAN 
expressions,  parentheses  can  be  used  to  override  the  usual  sequence  of 
evaluation.  Blanks  may  be  used  between  items  and  operator  symbols  to 
improve  readability  and  the  expression  may  extend  up  through  column  72  and 
continue  onto  the  next  card  with  no  continuation  marks.  The  end  of  the 
expression  must  be  indicated  by  a  dollar  sign  ($)  preceded  by  at  least  one 
blank. 
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H-Z-b-ii.  Calling  Sequences  of  WHELP  Matrix  Routines 

Because  it  may  on  occasion  be  desirable  either  to  call  a  WHELP 
subroutine  directly  or  to  know  its  calling  sequence  for  debugging  purposes, 
below  is  a  list  of  the  WHELP  operators  in  use  and  the  subroutine  calls  which 
will  result. 


Ope  ration 


Subroutine  Call 


C  =A+B  CA  LL  MATADD  (A,  N  RA  ,  NCA,  B,  NRB,  NCB,  C,  NR  DIMA, 

NRDIMB,  NRDIMC) 

C  =A  -  B  CALL  MATSUB  (A  ,  NR  A ,  NCA,  B,  NRB,  NCB,  C,  NRDIMA, 

NRDIMB,  NRDIMC) 


C=A*B 

T 

C=Ai*B 


C=Scalar*B 


CA  LL  MATMAT(A,  NRA,  NCA,  B,  NRB,  NCB,  C,  NRDIMA, 
NRDIMB,  NRDIMC) 

CALL  TRNSML(A,  NRA,  NCA,  B,  NRB,  NCB,  C,  NRDIMA, 
NRDIMB,  NRDIMC) 

CA  LL  SCAMAT  (SCALAR,  B,  NRB,  NCB,  C,  NRDIMB, 
NRDIMC) 


C  =  B/Scalar 
T 

C=A  1 
B  =  0. 

C=A.  B 

C=A  x  B 

B=A 

B^-A 

A=INVERSE  (A ) 
A  =IDENT  (A ) 


CALL  SCAMAT (1,  0/SCALAR,  B,  NRB,  NCB,  C,  NRDIMB, 
NRDIMC) 

CALL  T  R  NS  PS  (A ,  NRA,  NCA,  C) 

CALL  MATZRG5  (B,  NRB,  NCB,  NRDIMB) 

CA  LL  TRNSML(A,  NRA,  NCA,  B,  NRB,  NCB,  C,  NRDIMA, 
NRDIMB,  NRDIMC) 

CALL  C  R0SS  (A  ,  B,  C) 

CALL  MOVE  (A,  NRA  -NCA,  B) 

CALL  NEGATE  (A,  NRA*NCA,  B) 

CA  LL  MA TINV (A,  NRA,  0,  0,  DET,  NRDIMA) 

CALL  IDENT  (NRA  ,  A ) 


H  -  I  1 


where: 


A,  B,  and  C  are  WHELP  arrays  of  conformable  size  for  the  opera¬ 
tions  indicated. 

NRA  and  NRB  are  the  number  of  rows  of  A  and  B  being  used. 

NCA  and  NCB  are  the  number  of  columns  of  A  and  B  being  used. 

NR  DIMA,  NRDIMB,  and  NRDIMC  are  the  fixed  number  of  rows  for 
which  A,  B,  and  C  are  dimensioned. 

DET  is  the  determinant  of  matrix  A. 

H-2-c.  Special  Features 

•  A  =  0.  $  This  zeros  out  A  where  A  is  any  WHELP  variable. 

•  The  remainder  of  a  card  following  a  dollar  sign  may  be  used 
for  comment.  For  example: 

A  =  B  +  ERROR  $  ADD  ERROR  VECTOR 


•  WHELP  variables  may  be  set  equal  to  strings  of  FORTRAN 
expressions.  The  format  is 


variable  =  $  element,  $  element,,  $.  .  .  $  Element  $$ 

1  c  n 


where 


e  lement. 

i 


a  constant 
or 

any  legal  FORTRAN  expression 
or 

*k  (where  k  is  an  integer  constant  denoting 
how  many  times  the  expression  following  is 
to  be  repeated) 
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EXAMPLE: 

-^DECLARE  M(2,2)  $ 


M  =  $  SIN (T)  $  CO)S(T)  $  *2  $  ALPHA  +  R  $$ 
produces  the  result: 

M(l,  1)  =  SIN  (T ) 

M(2,  1)  =  C0S(T) 

M(l,  2)  =  ALPHA  +  R 
M(2,  2)  =  ALPHA  +  R 


Note : 

1.  Element  may  not  be  a  WHELP  expression. 

2.  Matrix  variables  are  stored  Ly  column  and  therefore  must  be 
listed  by  column. 

3.  Like  the  FORTRAN  DATA  statement,  no  elements  may  be 
skipped  and  the  number  of  elements  must  not  exceed  the  total 
size  of  the  variable. 


4.  When  used  with  variable  dimension  WHELP  (see  below),  data 
will  be  packed  into  the  first  N*M  elements  of  an  array. 

H-3.  VARIABLE  DIMENSION  WHELP 

WHELP  may  also  be  used  with  matrices  having  variable  dimen¬ 
sions:  the  WHELP  equations  are  written  in  exactly  the  same  manner  as  they 
are  for  fixed  dimension  WHELP,  and  special  forms  of  *SAMESIZE,  ^IDECLARE 
and  *INFQ)RM  are  used  to  declare  the  variable  size  matrices.  These  special 
forms  will  be  explained  below,  but  since  successful  use  of  variable  dimension 
WHELP  depends  upon  an  understanding  of  how  WHELP  stores  and  transmits 
data  for  matrices  of  variable  dimensions,  this  will  be  discussed  first.  It  is 
strongly  urged  that  the  user  carefully  observe  the  constraints  on  data  storage 
imposed  and  implied  for  variable  dimension  WHELP  and  also  that  thorough 
printout  and  testing  be  done  during  program  development. 
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H-3-a.  Data  Storage  and  Transmission 

We  are  accustomed  to  thinking  of  matrices  in  FORTRAN  as  having 
several  dimensions,  but  F0RTRAN  does  not  actually  store  a  matrix  in  a  two- 
dimensional  "slot":  It  stores  it,  column  by  column,  in  a  continuous  string. 
Thus  a  simple  two-dimensional  matrix  MAT(3,  3),  which  we  represent 
mathematically  as 


'I.  4.  7.- 

2.  5.  8. 

.3.  6.  9._ 

is  in  fact  stored  like  this: 

MAT  (1)  =  1. 

MAT  (2)  =  2. 

MAT  (3)  =  3. 

MAT  (4)  =  4. 

MAT  (5)  =  5. 

MAT{6)  =  6. 

MAT  (7)  =  7. 

MAT  (8)  =  8. 

MAT  (9)  =  9. 

As  long  as  full  matrices  are  used  with  WHELP  (or  F0RTRAN),  the  only 
commonly  encountered  implication  of  this  is  in  the  use  of  data  statements  to 
set  matrix  elements,  where  one  must  remember  to  list  data  by  columns  rather 
than  by  rows. 

Furthermore,  when  we  wish  to  deal  with  some  variable  size  subset  of  a 
matrix,  for  example,  if  we  want  to  use  the  above  MAT(3,  3)  as  MAT(N,  M) 
where  N  =  2  and  M  =  2,  then  we  are  accustomed  to  thinking  of  our  data  storage 
like  this : 
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MAT  (N,  M) 

'1. 

4.  I 

7.  ' 

2. 

1 

.i-j 

8. 

.3. 

6. 

9.  . 

MAT  (3,3) 


or  MAT (N,  M)  =  MAT(l.l)  =  1. 

MAT  (2,  1)  =  2. 
MAT (1,  2)  =  4. 
MAT (2,  2)  =  5. 


where  MAT  (N,  M)  occupies  the  first  2x2  positions  in  MAT (3,  3). 

WHELP,  however,  assumes  that  the  data  in  MAT  (N,  M)  is  stored 
in  the  first  N  x  M  locations  of  MAT(3,  3)  as  follows: 


MAT(N,  M) 


fl.  4.  1 

7.* 

or  MAT  (N,  M)  =  MAT(l) 

=  1. 

2  r7 

2.  j  0. 

1 

8. 

MAT  (2) 

=  2. 

.3.  j  6. 

9.. 

MAT  (3) 

=  3. 

MAT  (3,3) 

MAT(4) 

=  4. 

In  other  words,  WHELP  always  assumes  that  the  N  x  M  elements 
of  a  variable  dimension  matrix  are  stored  "packed",  one  immediately  after 
the  other,  by  columns,  in  the  storage  space  allotted  for  the  full  maximum  size 
of  the  array.  Whether  it  is  operating  on  a  matrix  or  storing  the  results  of  an 
operation  into  a  matrix,  it  will  use  the  first  N  x  M  elements,  not  the  first  N 
rows  and  M  columns. 

Therefore,  the  user  must  always  be  certain  that  arrays  to  be 
operated  on  are  stored  "packed"  and  that  if  WHELP  arrays  are  to  be  printed 
or  otherwise  used  in  FORTRAN  format,  they  must  be  "unpacked"  by  the  print 
statement  or  some  other  means.  Two  subroutines  are  included  in  the  WHELP 
library  to  aid  the  user  in  changing  from  FORTRAN  matrix  format  ("unpacked") 
to  WHELP  matrix  format  ("packed")  and  vice  versa.  They  are  explained 
below  in  Appendix  H-3-c. 


H-  15 


H-3-b.  Declaring  Variable  Dimension  WHELP  Variables 

Special  forms  of  *SAMESIZE,  ^IDECLARE,  and  *INF0RM 
accomplish  the  task  of  activating  variable  dimension  WHELP.  The  formats 
are  the  same  as  for  fixed  dimension  WHELP  except  that  n  and  m  (row  and 
column  dimensions)  can  take  either  of  two  forms: 

{integer  name /integer  constant 
or 

integer  constant 


where 

integer  name  must  be  1-3  characters,  beginning  with  a  letter 
integer  constant  is  the  maximum  size 


EXAMPLES: 

*SAMESIZE 

N/10 

M/20  A  B  $ 

*SAMESIZE 

20 

M/10  C  D  E  $ 

*IDECLARE 

A (N/ 10,  M/20)  B(10,  M/20)  C(5,  5)  $ 

*IDE  CLARE 

Z 

X(K/30)  Y  (K/20)  $ 

-INFORM 

N/20 

M/20  FM(1,1,K)  $ 

*INF<2RM 

L/20 

PT (1,  J)  $ 

The  following  statements  apply  to  all  three  declaration  forms: 

1.  FORTRAN  subroutine  calls  generated  by  WHELP  will  always 
use  the  letters  (if  any)  and  result  in  variable -dimension 
computations,  assuming  data  to  be  used  is  packed  and 
producing  packed  results. 

2.  The  numbers  given  as  dimensions  determine  the  maximum 
size  of  the  arrays  declared. 

3.  WHELP  checks  to  see  that  the  maximum  dimensions  of 
arrays  are  conformable  for  the  operations  indicated  in  an 
expression,  but  it  does  not  check  the  variable  dimensions 
to  ensure  they  are  less  than  the  maximum  dimensions  nor 
does  it  check  to  ensure  that  they  result  in  conformable 
matrices . 
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Beyond  this  lie  some  important  differences  in  how  the  three 
statements  may  be  used,  due  to  the  following  facts: 

1.  FORTRAN  requires  that  the  maximum  size  of  an  array 
must  be  stated  before  any  subset  of  the  array  may  be 
referenced. 

2.  If  an  array  is  to  have  variable  dimensions  within  a  sub¬ 
routine,  then  the  integer  variable  names  representing 
those  dimensions,  as  well  as  the  array  name,  must  be 
part  of  the  argument  list  of  the  subroutine. 

3.  *SAMESIZE,  *IDEC  LA  RE,  and  *INF0RM  are  all  trans  lated 
differently  by  the  WHELP  precompiler: 


*SAMESIZE  N/10  M/20  A  B  $ 

produces 

REAL  A (10,  20),  B(10,  20) 


whereas 


*IDE  CLARE  A  (N/10,  N/20)  B(10,M/20)  $ 

produces 

REAL  A(N,M),  B(10,M) 


and 


*INF0RM  N/20  PS(1,J)  $ 

produces 

(no  declaration  statement) 


Based  on  these  differences,  some  general  (though  by  no  means 
comprehensive)  guidelines  for  use  of  *SAMESIZE,  *IDECLARE,  and 
*INF0RM  may  be  suggested: 
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... 


1. 


*IDECLARE  may  not  be  used  in  a  main  program.  (Use 
*SAMESIZE.  ) 


2.  If  the  array  name  and  its  variable  dimensions  are  not  among 
the  subroutines  arguments,  *IDECLARE  may  not  be  used  to 
declare  the  array  in  a  subroutine. 

3.  *SAMESIZE  may  be  used  in  any  routine. 

4.  If  an  array  has  already  been  dimensioned  within  a  routine  by 
any  means,  ^INFORM  may  be  used  to  declare  the  entire  array 
or  any  totally  contiguous  subset  of  it  as  a  WHELP  variable. 
Remember,  though,  that  the  data  may  have  to  be  packed  if  it 
has  been  stored  in  FORTRAN  format. 

5.  Use  of  *INF0RM  to  declare  variable  dimensioned  subsets  of 
arrays  is  extremely  error  prone  due  to  the  contiguity  constraint 
and  should  be  used  only  with  great  care. 

6.  .  Results  should  be  checked  carefully,,  preferably  on  simple 

test  data,  as  many  possible  errors  will  not  produce  any 
warnings,  just  bad  results. 


EXAMPLE: 

PROGRAM  TEST WH (INPUT,  OUTPUT,  TAPE5=INPUT,  TAPE6=0UTPUT) 

*SAMESIZE  N/5  M/5  A  B  C  $ 

N  =  3 
M  =  3 

A  =  0.  $  ZEROES  0UT  FIRST  N  *  M  ELEMENTS  0F  A. 

B  =  IDENT(B)  $  CREATES  N  ORDER  IDENT  MATRIX,  PACKED. 

C  STORE  DATA  INTO  A  IN  PACKED  FORMAT 
A  =  $  *3  $  1.  $  *3  $  2.  $  *3  $  3.  $$ 

C  =  A  +  B  $ 

CALL  VARDIM(A,  B,C,N,  M) 

CALL  VARDIM2 (A ,  B,  C,  N,  M) 

END 
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SUBROUTINE  VARDIM(A2,  B2,  C2,  N,  M) 
*SAMESIZE  N/5  M/5  A2  B2  C2  $ 


C2  =  A2  +  B2  $ 

RETURN 

END 

SUBROUTINE  VARDIM2(A3,  B3,  C3,  N,  M) 

*IDEC  LARE  A3(N/5,  M/5)  B3(N/5,M/5)  C3(N/5,M/5)  $ 

C3  =  A3  +  B3  $ 

RETURN 

END 

In  this  example  C,  C2,  and  C3  will  all  wind  up  with  the  same  result, 
packed  into  the  first  N  *  M  locations.  Note  that  SUBROUTINE  VARDIM 
uses  *SAMESIZE  and  SUBROUTINE  VARDIM2  uses  *IDECLARE,  but 
results  are  identical. 


H-3-c.  Packing  and  Unpacking  Arrays 

To  aid  the  user  in  changing  from  FORTRAN  matrix  format 
("unpacked")  to  WHELP  matrix  format  ("packed")  and  vice  versa,  two 
subroutines  are  included  in  the  WHELP  library.  They  are  called  by: 


CALL  FT0WLP(A,  NR0WS,  NC0LS,  NDIMA) 

and 

CALL  WLPTOF(A,  NR0  WS,  NCOLS,  NDIMA) 


where 

A  is  the  matrix  to  be  packed  (unpacked) 

NR0WS  is  the  variable  row  dimension 

NCOLS  is  the  variable  column  dimension 

NDIMA  is  the  maximum  number  of  rows  for  which  A  is 
dimensioned 
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EXAMPLES: 

1.  *SAMESIZE  N/5  3  C  D  E  $ 

C  DEFINE  ELEMENTS  0 F  C  IN  FORTRAN  FORMAT 


C(1, 1) 

=  1. 

$ 

C(l,  2)  =  4. 

$ 

C(l,  3) 

C(2,  1) 

=  2. 

$ 

C(2,2)  =  5. 

$ 

C  (2,  3) 

C  (3,  1) 

=  3. 

$ 

C(3,  2)  =  6. 

$ 

C(3,  3) 

C  (4,  1) 

=  0. 

$ 

C (4,  2)  =  0. 

$ 

C  (4,  3) 

C(5,  1) 

N  =  3 

=  0. 

$ 

C(5, 2)  =  0. 

$ 

C (5,  3) 

C  PUT  C  INTO  WHELP  FORMAT 
CALL  FT0WLP(C,  N,  3,  5) 

C  SET  E  =  C (N,  3)  X  AN  N  X  3  IDENTITY  MATRIX 
E  =  C  *  IDE  NT  (D)  $ 

Notes: 

a.  In  this  example  the  elements  of  C  are  all  defined  in  FORTRAN 
format,  so  before  C(N,  3)  can  be  used  in  a  WHELP  equation, 
the  data  must  be  packed  into  the  first  N  x  3  elements  of  C. 
IDENT(D)  will  be  a  packed  matrix,  as  will  matrix  E  since  they 
are  the  results  of  WHELP  operations. 

b.  The  packing  operation  will  destroy  values  previously  stored  in 
C(4,  1)  through  C(4,  2).  A  subsequent  unpack  would  leave 
changed  values  in  C(4,  1),  C(5,  1)  and  C(4,  2). 

2.  SUBROUTINE  XYZ(A,  N,  M,  NDIMA) 

DIMENSION  A  (NDIMA,  1) 

*INF0RM  N/20  M/3  A  $ 

C  PUT  INTO  WHELP  FORMAT 

CALL  FTOWLP(A,N,M,  NDIMA) 

C  FORM  A  *  A -TRANSPOSE  FOR  A  BEING  N  X  M 
A  =  A  *  A,  $ 

C  RETURN  RESULT  IN  FORTRAN  FORMAT 
CALL  WLPT0F(A,  N,  M,  NDIMA) 

RETURN 

END 


H-20 


In  this  example  matrix  A  was  not  a  variable  dimension  WHELP 
array  in  the  calling  program,  so  it  comes  to  the  subroutine  un¬ 
packed  and  returns  unpacked,  but  must  be  packed  in  order  for 
the  variable  dimension  WHELP  equation  to  be  executed  correctly. 
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general  rules  .  3-2 

numbe  rof . 7-1,  7-2 

with  WHELP .  3-3,  3-4 

examples .  3-4 

DERIVS  (subroutine) . 3-1,  3-2 

calls  to  during  switch  processing . 5-16,  5-17 

example . 2-8 

format . D-6 

printing  from . 6-12 

when  written  .  D-5 

*DERIVS .  3-1 

DIFTAB  (subroutine)  .  D-9 

dimensions  (maximum)  . .  F-4 

discontinuities . . . .  4-8,  5-1  to  5-19 

accuracy  control  .  5-18 

debugging . .  .  G-2 

detecting  timing .  5-16 

printing  at . 5-17,  5-18 

sequence  of  events .  5-15  to  5-17 

several  in  series  . 5-17 

See  also  *SWTCH,  *SWMEM,  *EVENT  . 
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discontinuous  driving  function  (see  *SWTCH) 

discrete  inputs . 4-7,  4-8 

DY  .  3-1  to  3-4,  C-l,  C-4 


*ENDDERIVS . 

*ENDIC  (see  *ICC(JMP) 
*END(?UT  (see  OUTPUT) 
EPS . 


*EPS .  4-3 

example . 4-4 

ESPCT  L  (subroutine)  . 4-2,  D-8 

flowchart  . . . . . D-13  to  D-15 

ESPII  (subroutine) .  D-8,  D-10 

flowchart  . D-12 

ESPlijT  (subroutine)  . .  D-8 

ESPPC  (subroutine) 

equations  . E-3,  E-4 

flowchart  . D-20  to  D-22 

See  also  integration  package.  Predictor  Corrector. 

ESPPLOT  (program)  .  6-13  to  6-16 

example . 6-15,  6-16 

ESPRINT  (subroutine) . D-8 

ESPRK2  (subroutine) 

flowchart  . D-18,  D-19 

equations  . E-l,  E-2 

See  also  integration  package,  Runge-Kutta. 

ESPRK4  (subroutine) 

flowchart  . D-18,  D-19 

equations  . E-2,  E-3 

See  also  integration  package,  Runge-Kutta. 

EVENTS  (subroutine). . 5-10,  5-11,  5-16 


FILBUF  (subroutine)  . 

files 

file  usage  for  ESP  without  WHELP . 

file  usage  for  ESP  with  WHELP  . 

See  also  TAPE11,  TAPE12,  TAPE15,  datafiles. 

FIRST  P . 

See  also  ESPRK4,  flowchart. 


4-7,  4-8,  C-6 
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FIX  ST  P . 

flowcharts  . 

ADAMS  . 

flowchart  conventions  .... 

ESPCTL  . 

ESPII  . 

ESP  package,  overview  .  .  . 

ESPPC  . 

ESPRK4  (ESPRK2) . 

formats,  card 

general  rules  . 

summary  . 

See  also  specific  card  names. 
FT(JWLP  (subroutine) . 


4-3,  4-6,  4-8,  C-7 
...  D-10  to  D-22 
.  .  .  .  D- 16,  D- 17 

. D-ll 

.  .  .  D- 13  to  D-15 

. D  -  x 

. D-10 

.  .  .  D-20  to  D-22 
.  .  .  D-  18,  D- 19 


.  .  A  -  1 
A-2,  A-3 


H- 19  to  H-21 


G 


GRAPH  (subroutine)  (see  *GRAPH) 


♦GRAPH  . 6-5  to  6-11 

examples . 2-4,  6-10 

use  in  multiple  runs . F-l  to  F-3 

graphic  output  (see  plotting) 

GRIDggg  .  6-8 


H 


H . 

See  also  stepsize. 

Hamming  (see  integration  package,  predictor-corrector) 

HEAD . 

♦HEADINGS  . 

HHMAX . 

HHMIN . 

HMAX  . 

HMAXMN  (common  block)  . . 

HMIN . 

HP  . 

hprinti  . 

See  also  HP. 

HSW . 

*HSW . 

examples  . 


4-2  to  4-6,  C-7 


.  C-10 

.  6-2,  7-8 

.  C-4 

.  C-4 

.  4-3,  C-7 

.  C-4 

4-3,  4-10,  C-7,  G-2 

.  4-3,  C-7 

.  7-1,  7-2 

.  C-5 

.  5-16,  5-18 

.  %JL3 
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HSWE  . 

*HSWE . 

examples  . 

HSWM  . 

*HSWM . 

examples  . 

hysteresis  (see  *SWMEM) 


I 


ICC(jMP  (subroutine) . 

example  of  dummy  .  . . 

format . 

printing  from . . 

when  written . 

*ICC0MP . 

example . 

*IDEC  L*ARE 

fixed  dimension . 

example . 

variable  dimension . 

examples . 

IDENT  (subroutine) . 

IE  VENT . 

IF0RM . 

Ill-conditioned  system  ........ 

*IMAX  . 

IMS0RC  (file)  . 

independent  variable  (see  T) 

*INF0RM 

fixed  dimension . 

examples . .  .  .  .  . 

variable  dimension . 

examples . 

initial  conditions  (see  *IV,  *ICC0MP) 

inputs  . . 

initial  conditions  . 

miscellaneous . 

run  times  . . 

to  control  accuracy . 

user  cards  . 

user  files  . . 

user  parameters  . 


.  .  C-5 

5-16,  5-18 
.  .  5-19 

.  .  C-5 

5-16,  5-18 
.  .  5-19 


.  7-5 

.  2-12 

......  D-7 

.  6-11 

.  D-5 

.  .  .  7-5,  7-6 

.  7-6 

.  .  .  H-3,  H-4 
H-4 

.  H-16  to  H-19 
.  .  H-16,  H- 19 

. H-ll 

.  .  .  5-11,  C-l 

. C-10 

4-9,  4-10,  G-2 

.  6-14 

.  .  .  D-2,  D-4 


.  .  .  .  H-5  to  H-7 

.  H-6 

.  .  .  H-16  to  H-19 
.  .  .  .  H-16,  H-20 

.  .  .  .  7-1  to  7-9 

7-2,  7-3,  7-5,  7-6 

.  7-7 

.  7-1,  7-2 

.  7-7 

.  7-6,  F-2 

. 7-6,  F-2 

........  7-3 
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integration  package . { 

Adams  integration . 

alternatives . . 

error  control  . 

Predictor-Corrector . 

Runge-Kutta . . . 

solution  accuracy . 

stepsize  control  .  .  . 

integrator  equations . 

INTERP  (see  INRKPC,  ADMNTP) 

interpolation  of  printed  values . .  . 

ISWTYP . 

*IV . 


. 4-1  to  4-10,  E-l  to  E-4 

.  .  .  .  4-4,  4-5,  E-l,  D- 16,  D-17 

.  4-1 

.  .  .  4-4  to  4-7,  4-10,  E-l  to  E-4 

4-9,  4-10,  E-3,  E-4,  D-ZO  to  D-22 

4-6  to  4-8,  E-l  to  E-3,  D-18,  D-19 

.  4-3 

.  4-2,  4-3 

. E-l  to  E-4 

.  6-4,  6-5 

.  C-8 

.  7-2,  7-3 


J 


JUNE  .  C-10 

job  control  cards  (see  control  cards) 

J START .  C-7 


K 


KC0UNT . 4-9,  D-21  to  D-22 

KFLAG .  C-7 


L 


large  simulations . F-4,  F-5 

LINES  .  C-10 


M 


MAIN  (PROGRAM) 

examples  .  2-7,  6-17 

format . D-6 

when  written .  D-4 

MATADD  (subroutine) . H-ll 

MATINV  (subroutine)  . H-ll 

MATMAT  (subroutine) . H-ll 

matrix  equations  (see  WHELP) 

MATSUB  (subroutine)  .  . . H-ll 

MATZRO  (subroutine) .  H-ll 
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MAX .  C-9 

MAXCHR .  C-6 

MAXCOL .  C-6 

MAXDER .  C-7 

MAXMEM .  C-9 

*MAXPL<jTS .  6-6,  6-7,  6-14 

MAXSWS .  C-9 

*METHtf>D . 4-1,  D-4 

MF .  C-7 

MISCEL  (common  block) .  C-5 

M0VE  (subroutine) . 3-3,  3-4,  H- 11 

multiple  runs . F-ltoF-3 

MXL .  C-10 


N 


NALTER . 

NCHNG  . 

NDISPR . 

NDISPR  (common  block) . 

NEGATE  (subroutine) . 

NEQ  (neg) . 

NEVE NT . 

*NEVENT  . 

NFIRST . 

NHEAD  . 

NL0CAL  . 

*NL0CAL . 

noise  inputs . 

nonlinearities  (see  discontinuities) 

N0PL0T . 

NOTIFY  (subroutine) . 

NPAGE  . 

NP0INT . 

NTAP11 . 

NUMSTP . 


.  C-8 

C  -  8 

5-16,  *5*- 17,’  5-18,  C-6 

. 5-18,  C-6 

. H-ll 

.  .  .  .  3-2,  7-1,  C-5 

.  C-9 

.  5-11 

.  C-6 

.  C-10 

.  C-9 

.  6-14 

.  4-8 

.  C-9 

. 5-11,  5-16 

. C-10 

.  C-9 

.  C-4 

.  C-4 


O 

0DES0L  (see  ESPPC,  ESPRK2,  ESPRK4,  ADAMS) 

0UT .  C-10 
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output  . . 

accuracy  ....  . 

alternatives . 

data  file . 

graphic . 

listing . 

initial  printout  example  .  . 
printer  plot  examples  .  .  . 
user  printout  examples  .  . 
plotted  (see  graphic) 

printed  . 

tape  (see  data  file) 

See  also  *PRINT,  *HEA DINGS. 

(OUTPUT  (subroutine) . . 

example . 

format . 

printing  from . 

when  written . 

♦Output . .  .  .  . 

examples . 

0VERLAY  . 

0VERLAY1 . 

overlays  (plotting) . 

from  multiple  runs  ...... 


. 6- 1  to  6- 17 

.  .  .  .  6-4,  6-5,  6-7 

.  6-1 

. 6-13  to  6-17 

. 6-5  to  6-11 

.  .  .  .  2- 16  to  2-  19 

.  2-16 

.  2-18,  2-19 

.  2-17 

6-1  to  6-5,  6-11,  6-12 


.  .  6-2 
.  .  2-11 
.  .  D-7 

.  .  6-12 
.  .  D-5 

6-3,  6-4 
6-4,  6-6 
.  .  6-8 
.  .  6-8 
.  .  6-8 
F-l  to  F-3 


P 


PACKER  (subroutine) . 

PAR . 

examples  . 

*PAR . 

example  s  . 

parameters  (see  *PAR) 

PARS  (common  block) . 

PD0T . 

plotting . 

accuracy  . 

alternatives . 

controlling  at  switch  times 

from  TAPE  11 . 

number  of  plots . 

overlays  from  different  runs 


.  D-9 

.  7- 3,  7-4,  C-6 

.  2-3,  7-4 

.  7-3,  7-4 

.  2-3,  7-4 

.  C-6 

.  .  6-6,  6-13,  C-l,  C-4 
6-5  to  6-11,  6-13  to  6-16 

.  6-7 

.  6-5 

. 5-17,  5-18 

. 6-13  to  6-16 

.  6-6 

. F-2,  F-3 
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storing  plot  data . 

using  *GRAPH . 

using  Program  ESPPLQiT 
PREC0MP  (Program)  .... 

cards  written  by . 

how  used  . 

sequence  of  events  .... 
termination  during  .... 

PRINT  . 

PRINT  (1) . 

examples . 

*PRINT . 

examples . 

printing . 

accuracy  . 

controlling  at  switch  times 

in  *ICC0MP . 

in  *DERIVS  . 

in  MAIN . 

in  *0UTPUT . 

user -formatted . 

print  interval  (see  HP,  hprinti) 
printout,  debugging . 


.  6-6 

. 6-5  to  6-11 

. 6-13  to  6-16 

.  .  .  7-6,  D-4  to  D-7 

. D-5  to  D-7 

. D-2,  D-3 

. D-4,  D-5 

.  G-l 

.  C-l 

. 6-5,  6-12 

.  6-12,  6-17 

.  .  .  .  6-1  to  6-4,  6-6 
.  .  2-4,  6-2,  6-3,  6-4 
6-1  to  6-5,  6-11,  6-12 

.  6-4,  6-5 

. 5-17,  5-18 

.  6-11 

.  6-12 

.  6-12 

.  6-12 

. 6-11,  6-12 


G-l,  G-2 


Q 


Q .  C-5 

*Q .  4-3 

example . 4-4 


R 


READIN  (common  block) . C-6 

READIT  (subroutine)  . D-4,  D-9 

reserved  names  . C-ll,  C-12 

*  RETURN .  7-9 

RKC0NT  (common  block) 

format  and  contents .  C-6 

use . 4-7 

*RUN  .  7-1,  7-2 

*RUNC . F-2,  F-3 

RUNGE-KUTTA  integration . . . 4-6  to  4-8 

run-time  routines . . . D-7  to  D-9 
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*SAMESIZE 

fixed  dimension . . 

example . . 

variable  dimension.  .... 

examples . . 

SCALE . 

SCAMAT  (subroutine)  . . 

SCR  I  (200) . 

SECNZR  (subroutine) . 

SIZExxyy . 

SKPFIL  (subroutine) . 

SMALL  . 

stacked  runs  . 

stepsize 

control 

by  Adams . 

by  Predictor-Corrector 
by  Runge-Kutta  .... 

by  user  . 

initial  selection . 

relationship  to  print  interval 
STFPAR  (common  block)  .... 

ST0P . 

*ST0P  . 

STPC0N  (common  block)  .... 
format  and  contents  .... 

use . 

subroutines  (see  specific  names) 

SWCHi  . 

examples . 

SWDBUG  . 

SWHPAR  (common  block)  .  .  .  . 

SWINIT  (subroutine) . 

SWINPT  (subroutine) . 

examples  . 

format . 

user-written  ........ 

^SWITCHES . 

SWMi . 

SWMEM . 

SWMEM(i,  4) . 


. H-3,  H 

.  H 

. H- 16  to  H- 

H-  16,  H- 18,  H-  19,  H- 

.  6 

.  H- 

. .  C 

.  D 

.  6 

.  D 

.  6 

. F-l  to  F 


4-4, 

4-9, 


4-2, 


6-4, 


C-l, 


4-2, 


.  .  5-3,  5-14,  5- 

.  5 

.  C 

. 5-9,  C 

.  D 

.  5-2,  5-12  to  5- 

.  2 

.  D 

.  5- 

.  5- 

.  .  .5-7,  5-14,  5- 

.  C 

5-7,  5-14,  5-15,  5- 


i 
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*SWMEM .  5-5  to  5-10,  5-12  to  5-18 

defining  inputs . 5-6,  5-12  to  5-14 

defining  output .  5-7,  5-14,  5-15 

defining  characteristics  . 5-7,  5-8,  5-9 

how  it  works .  5-15  to  5-18 

*SWMEMCNT .  5-14 

*SWMEMDATA . 5-6,  5-7  to  5-9 

SWMEMN  (subroutine) . 5-6,  5-12  to  5-14 

example  of  dummy .  2-10 

format .  D-7 

user-written . 5-14 

*SWMEMSET . 5-6,  5-10 

SWSET .  C-8 

SWTCH(i) . 5-3,  5-14,  5-15,  5-17,  C-8 

examples . 5-4 

*SWTCH  . 5-2  to  5-5,  5-12  to  5-18 

defining  inputs . 5-2  to  5-4,  5-12  to  5-14 

defining  output .  5-3,  5-4,  5-14,  5-15 

examples  . 2-3,  5-3,  5-4 

how  it  works .  5-15  to  5-18 

SWTCHE  (subroutine)  .  D-9 

SWTCHS  (common  block) . C-8,  C-9 


T . 3-2,  4-4,  4-6,  4-7,  4-9,  4-10,  C-l 

TO . 4-2,  C-3 

TAPE11 . 6-7,  6-13,  B-4,  B-5,  D-2,  D-3 

TAPE  12 . 7-6,  B-4,  B-5,  D-2,  D-3,  D-5,  F-2 

TAPE  15 . D-2,  D-3 

TERMCH .  C-6 

TF .  C-3 

def.,  length,  and  block  . .  C-3 

*TICKPLQ>TS .  6-14 

time,  start  and  stop . 7-1,  7-2 

TIMEIN  (subroutine) . D-9 

TIME<JUT  (subroutine) .  D-9 

TITLE  . C-10 

*TITLE .  7-8 

example . 2-5 

TQJDAY .  c-10 

TP  .  C-3 
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TRNSMUL  (subroutine)  . . . H-ll 

TRNSPS  (subroutine)  . .  H-ll 

TYPE  tl  (t2)(t3)  . 6-9,  6-10 


UNIP1  (common  block) .  C-9 

UNIP2  (common  block) .  C-10 


VALEVS . C-5 

VA  LMEM .  C-5 

VALUES  . C-l,  C-5 

in  EVENTS .  5-11 

in  SWINPT . 5-2,  5-14 

in  SWMEMN . 5-6,  5-14 

variables,  alphabetical  list  .  C-l 

VPL0T .  C-9 


WHELP . H-l  to  H-21 

data  storage . H-15,  H-15 

declaring  variables 

fixed  dimensions . H-2  to  H-7 

variable  dimensions  . . . H-16  to  H-19 

defining  derivatives  in  .  . .  3-3,  3-4 

how  it  works  .  .  .  . . H-l,  H-2 

matrix  routines  . . H-ll 

operator  symbols . H-10 

packing  data . H-19  to  H-21 

special  features  .  H-12 

writing  WHELP  expressions . .  . .  H-8  to  H-13 

See  also  *IDECLARE,  *SAMESIZE,  *  INFORM  . 

Wilkinson's  method  .  5-16 

WLPTOF  (subroutine)  .  . . H-19  to  H-21 


3-1  to  3-7, 


See  also  *IV. 
YPRNT  (100)  .  . 
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